diff --git a/.nojekyll b/.nojekyll index 2577ce2..a792f62 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -4268e2e2 \ No newline at end of file +1f3656a0 \ No newline at end of file diff --git a/01-introduccion.html b/01-introduccion.html index 65942b3..4b5efd7 100644 --- a/01-introduccion.html +++ b/01-introduccion.html @@ -287,18 +287,18 @@

Ejemp -B -chicos +A +grandes mejora B -chicos +grandes mejora -B -chicos +A +grandes mejora @@ -307,12 +307,12 @@

Ejemp mejora -B -chicos +A +grandes sin_mejora -B +A grandes mejora @@ -322,19 +322,19 @@

Ejemp mejora -A -grandes +B +chicos mejora A grandes -sin_mejora +mejora A grandes -sin_mejora +mejora @@ -578,8 +578,8 @@

Ejemp ", width = 200, height = 50)
-
- +
+

Es decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).

@@ -696,8 +696,8 @@

Eje ", width = 200, height = 50)
-
- +
+

Nótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas.

diff --git a/02-flujo-basico-2.html b/02-flujo-basico-2.html index 31220ef..31bcdbe 100644 --- a/02-flujo-basico-2.html +++ b/02-flujo-basico-2.html @@ -203,7 +203,7 @@

3 

3.1 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:

@@ -226,7 +226,7 @@

node [shape=plaintext] N Npos [label = <N<SUB>+</SUB>>] - Nobs + Nobs [label = <N<SUB>obs</SUB>>] #Nneg [label = <N<SUB>-</SUB>>] #sens #esp @@ -247,8 +247,8 @@

")#, 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.

@@ -290,7 +290,7 @@

3.1.3 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

+

Si llamamos a esta cantidad \(\theta_{obs}\), de forma que dada una muestra 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 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_{-}}\]

@@ -374,7 +374,12 @@

geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + facet_wrap(~rep)
-

+
+
+

+
Figura 3.1: Verificación a priori
+
+

Contrasta con lo que pasaría si usaramos el modelo sin considerar fuentes de error:

@@ -402,7 +407,12 @@

geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + facet_wrap(~rep)
-

+
+
+

+
Figura 3.2: Verificación a priori fallida (modelo incorrecto)
+
+

Este resultado está lejos de ser aceptable.

@@ -465,8 +475,8 @@

")#, 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:

@@ -503,7 +513,7 @@ 

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"))
resumen |> select(variable, mean, q5, q95)
@@ -571,12 +581,12 @@

-
ggplot(sims, aes(x = p)) + 
+
ggplot(sims, aes(x = theta)) + 
   geom_histogram()

 

 
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

@@ -588,7 +598,7 @@ 

-
ggplot(sims, aes(x = esp, y = p)) + geom_point() +
+
ggplot(sims, aes(x = esp, y = theta)) + geom_point() +
   xlab("Especificidad del kit") + ylab("Prevalencia") + geom_smooth()

 

 
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

@@ -599,7 +609,7 @@ 
\(p\) es considerable: va desde valores cercanos a 0 hasta valores alrededor de 0.025-0.03.

+
Y notamos que aún con una muestra relativamente grande, el rango de \(\theta\) es considerable: va desde valores cercanos a 0 hasta valores alrededor de 0.025-0.03.

 
 
 
diff --git a/02-flujo-basico-2_files/figure-html/unnamed-chunk-8-1.png b/02-flujo-basico-2_files/figure-html/fig-apriori-check-1.png
similarity index 100%
rename from 02-flujo-basico-2_files/figure-html/unnamed-chunk-8-1.png
rename to 02-flujo-basico-2_files/figure-html/fig-apriori-check-1.png
diff --git a/02-flujo-basico-2_files/figure-html/unnamed-chunk-10-1.png b/02-flujo-basico-2_files/figure-html/fig-apriori-falla-1.png
similarity index 100%
rename from 02-flujo-basico-2_files/figure-html/unnamed-chunk-10-1.png
rename to 02-flujo-basico-2_files/figure-html/fig-apriori-falla-1.png
diff --git a/02-flujo-basico-2_files/figure-html/unnamed-chunk-17-1.png b/02-flujo-basico-2_files/figure-html/unnamed-chunk-17-1.png
index d223478..f615f5d 100644
Binary files a/02-flujo-basico-2_files/figure-html/unnamed-chunk-17-1.png and b/02-flujo-basico-2_files/figure-html/unnamed-chunk-17-1.png differ
diff --git a/02-flujo-basico.html b/02-flujo-basico.html
index 23cdf5d..3060838 100644
--- a/02-flujo-basico.html
+++ b/02-flujo-basico.html
@@ -263,8 +263,8 @@ 
", width = 300, height = 100)

 
 

-

-
+

+
 

 

 
Que también podríamos simplificar (suponiendo la \(N\) fija y conocida, pues \(N_+\) y \(M\) dan \(N_{-}\)) como:

@@ -297,8 +297,8 @@ 
", 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.

@@ -989,7 +989,7 @@

@@ -1288,7 +1288,10 @@

<

\[\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.

+

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 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. diff --git a/search.json b/search.json index 9d9871b..8d08fe2 100644 --- a/search.json +++ b/search.json @@ -11,7 +11,7 @@ "href": "01-introduccion.html#diagramas-causales", "title": "1  Introducción", "section": "1.1 Diagramas causales", - "text": "1.1 Diagramas causales\nEn primer lugar, observamos (McElreath (2020)):\n\n\n\n\n\n\nCausas y mecanismos\n\n\n\nLas razones de cómo hacemos análisis estadístico (que procedimiento o algoritmo seleccionamos, por ejemplo) en un problema dado no están en los datos observados, las causas de los datos.\n\n\nLas causas de los datos no pueden extrarse de los datos solamente. Muchas veces nos referimos a las causas de los datos como el proceso generador de los datos: esto incluye aspectos del fenómeno que nos interesa (ciencia o proceso de negocios, etc.), así como el proceso de observación (muestras, valores no observados, etc.).\nConsideremos un ejemplo simple para ilustrar este primer principio:\n\nEjemplo (cálculos renales)\nEste es un estudio real acerca de tratamientos para cálculos renales (Julious y Mullee (1994)). Pacientes se asignaron de una forma no controlada a dos tipos de tratamientos para reducir cálculos renales. Para cada paciente, conocemos el tipo de ćalculos que tenía (grandes o chicos) y si el tratamiento tuvo éxito o no.\nLa tabla original tiene 700 renglones (cada renglón es un paciente)\n\ncalculos <- read_csv(\"../datos/kidney_stone_data.csv\")\nnames(calculos) <- c(\"tratamiento\", \"tamaño\", \"éxito\")\ncalculos <- calculos |> \n mutate(tamaño = ifelse(tamaño == \"large\", \"grandes\", \"chicos\")) |> \n mutate(resultado = ifelse(éxito == 1, \"mejora\", \"sin_mejora\")) |> \n select(tratamiento, tamaño, resultado)\nnrow(calculos)\n\n[1] 700\n\n\ny se ve como sigue (muestreamos algunos renglones):\n\ncalculos |> \n sample_n(10) |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\n\n\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\nsin_mejora\n\n\nB\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nA\ngrandes\nsin_mejora\n\n\n\n\n\n\n\nAunque estos datos contienen información de 700 pacientes, los datos pueden resumirse sin pérdida de información contando como sigue:\n\ncalculos_agregada <- calculos |> \n group_by(tratamiento, tamaño, resultado) |> \n count()\ncalculos_agregada |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\nn\n\n\n\n\nA\nchicos\nmejora\n81\n\n\nA\nchicos\nsin_mejora\n6\n\n\nA\ngrandes\nmejora\n192\n\n\nA\ngrandes\nsin_mejora\n71\n\n\nB\nchicos\nmejora\n234\n\n\nB\nchicos\nsin_mejora\n36\n\n\nB\ngrandes\nmejora\n55\n\n\nB\ngrandes\nsin_mejora\n25\n\n\n\n\n\n\n\nComo en este caso nos interesa principalmente la tasa de éxito de cada tratamiento, podemos mejorar mostrando como sigue:\n\ncalculos_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, tamaño, total, prop_mejora) |> \n arrange(tamaño) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\ntotal\nprop_mejora\n\n\n\n\nA\nchicos\n87\n0.93\n\n\nB\nchicos\n270\n0.87\n\n\nA\ngrandes\n263\n0.73\n\n\nB\ngrandes\n80\n0.69\n\n\n\n\n\n\n\nEsta tabla descriptiva es una reescritura de los datos, y no hemos resumido nada todavía. Pero es apropiada para empezar a contestar la pregunta:\n\n¿Qué indican estos datos acerca de qué tratamiento es mejor? ¿Acerca del tamaño de cálculos grandes o chicos?\n\nSupongamos que otro analista decide comparar los pacientes que recibieron cada tratamiento, ignorando la variable de tamaño:\n\ncalculos |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\ny parece ser que el tratamiento \\(B\\) es mejor que el \\(A\\). Esta es una paradoja (un ejemplo de la paradoja de Simpson) . Si un médico no sabe que tipo de cálculos tiene el paciente, ¿entonces debería recetar \\(B\\)? ¿Si sabe debería recetar \\(A\\)? Esta discusión parece no tener mucho sentido.\nPodemos investigar por qué está pasando esto considerando la siguiente tabla, que solo examina cómo se asignó el tratamiento dependiendo del tipo de cálculos de cada paciente:\n\ncalculos |> group_by(tratamiento, tamaño) |> count() |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nn\n\n\n\n\nA\nchicos\n87\n\n\nA\ngrandes\n263\n\n\nB\nchicos\n270\n\n\nB\ngrandes\n80\n\n\n\n\n\n\n\nNuestra hipótesis aquí es que la decisión de qué tratamiento usar depende del tamaño de los cálculos. En este caso, hay una decisión pues A es una cirugía y B es un procedimiento menos invasivo, y se prefiere utilizar el tratamiento \\(A\\) para cálculos grandes, y \\(B\\) para cálculos chicos. Esto quiere decir que en la tabla total el tratamiento \\(A\\) está en desventaja porque se usa en casos más difíciles, pero el tratamiento \\(A\\) parece ser en general mejor. La razón es probablemente un proceso de optimización de recursos y riesgo que hacen los doctores.\n\nEn este caso, una mejor respuesta a la pregunta de qué tratamiento es mejor es la que presenta los datos desagregados.\nLa tabla desagregada de asignación del tratamiento nos informa acerca de cómo se está distribuyendo el tratamiento en los pacientes.\n\n\n\n\n\n\n\nNota\n\n\n\nLos resúmenes descriptivos acompañados de hipótesis causales acerca del proceso generador de datos, nos guía hacia descripciones interpretables de los datos.\n\n\nLas explicaciones no son tan simples y, otra vez, interviene el comportamiento de doctores, tratamientos, y distintos tipos de padecimientos.\nPodemos codificar la información causal con un diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n T \n M \n C\n edge [minlen = 3]\n T -> M\n C -> T\n C -> M\n{ rank = same; M; T }\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).\n\n\nEjemplo (cálculos renales 2)\nContrastemos el ejemplo anterior usando exactamente la misma tabla de datos, pero con el supuesto de un proceso generador diferente. En este caso, los tratamientos son para mejorar alguna enfermedad del corazón. Sabemos que parte del efecto de este tratamiento ocurre gracias a una baja en presión arterial de los pacientes, así que después de administrar el tratamiento, se toma la presión arterial de los pacientes. Ahora tenemos la tabla agregada y desagregada como sigue:\n\ncorazon <- calculos |> \n select(tratamiento, presión = tamaño, resultado) |> \n mutate(presión = ifelse(presión == \"grandes\", \"alta\", \"baja\"))\ncorazon_agregada <- corazon |> \n group_by(tratamiento, presión, resultado) |> \n count()\ncorazon_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, presión, total, prop_mejora) |> \n arrange(presión) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\npresión\ntotal\nprop_mejora\n\n\n\n\nA\nalta\n263\n0.73\n\n\nB\nalta\n80\n0.69\n\n\nA\nbaja\n87\n0.93\n\n\nB\nbaja\n270\n0.87\n\n\n\n\n\n\n\n\ncorazon |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\n¿Cuál creemos que es el mejor tratamiento en este caso? ¿Deberíamos usar la tabla agregada o la desagregada por presión?\n\nEn este caso, la tabla agregada es más apropiada (B es mejor tratamiento).\nLa razón es que presión en este caso es una consecuencia de tomar el tratamiento, y como las tablas muestran, B es más exitoso en bajar la presión de los pacientes.\nSi sólo comparamos dentro de los grupos de presión baja o de presión alta, ignoramos lo más importante del tratamiento en la probabilidad de mejorar.\n\nNuestros supuestos causales podemos mostrarlos con el siguiente diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n P\n T \n M \n edge [minlen = 3]\n T -> P\n P -> M\n T -> M\n{ rank = same; M; T}\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nNótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas." + "text": "1.1 Diagramas causales\nEn primer lugar, observamos (McElreath (2020)):\n\n\n\n\n\n\nCausas y mecanismos\n\n\n\nLas razones de cómo hacemos análisis estadístico (que procedimiento o algoritmo seleccionamos, por ejemplo) en un problema dado no están en los datos observados, las causas de los datos.\n\n\nLas causas de los datos no pueden extrarse de los datos solamente. Muchas veces nos referimos a las causas de los datos como el proceso generador de los datos: esto incluye aspectos del fenómeno que nos interesa (ciencia o proceso de negocios, etc.), así como el proceso de observación (muestras, valores no observados, etc.).\nConsideremos un ejemplo simple para ilustrar este primer principio:\n\nEjemplo (cálculos renales)\nEste es un estudio real acerca de tratamientos para cálculos renales (Julious y Mullee (1994)). Pacientes se asignaron de una forma no controlada a dos tipos de tratamientos para reducir cálculos renales. Para cada paciente, conocemos el tipo de ćalculos que tenía (grandes o chicos) y si el tratamiento tuvo éxito o no.\nLa tabla original tiene 700 renglones (cada renglón es un paciente)\n\ncalculos <- read_csv(\"../datos/kidney_stone_data.csv\")\nnames(calculos) <- c(\"tratamiento\", \"tamaño\", \"éxito\")\ncalculos <- calculos |> \n mutate(tamaño = ifelse(tamaño == \"large\", \"grandes\", \"chicos\")) |> \n mutate(resultado = ifelse(éxito == 1, \"mejora\", \"sin_mejora\")) |> \n select(tratamiento, tamaño, resultado)\nnrow(calculos)\n\n[1] 700\n\n\ny se ve como sigue (muestreamos algunos renglones):\n\ncalculos |> \n sample_n(10) |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\n\n\n\n\nA\ngrandes\nmejora\n\n\nB\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\n\n\n\n\n\nAunque estos datos contienen información de 700 pacientes, los datos pueden resumirse sin pérdida de información contando como sigue:\n\ncalculos_agregada <- calculos |> \n group_by(tratamiento, tamaño, resultado) |> \n count()\ncalculos_agregada |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\nn\n\n\n\n\nA\nchicos\nmejora\n81\n\n\nA\nchicos\nsin_mejora\n6\n\n\nA\ngrandes\nmejora\n192\n\n\nA\ngrandes\nsin_mejora\n71\n\n\nB\nchicos\nmejora\n234\n\n\nB\nchicos\nsin_mejora\n36\n\n\nB\ngrandes\nmejora\n55\n\n\nB\ngrandes\nsin_mejora\n25\n\n\n\n\n\n\n\nComo en este caso nos interesa principalmente la tasa de éxito de cada tratamiento, podemos mejorar mostrando como sigue:\n\ncalculos_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, tamaño, total, prop_mejora) |> \n arrange(tamaño) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\ntotal\nprop_mejora\n\n\n\n\nA\nchicos\n87\n0.93\n\n\nB\nchicos\n270\n0.87\n\n\nA\ngrandes\n263\n0.73\n\n\nB\ngrandes\n80\n0.69\n\n\n\n\n\n\n\nEsta tabla descriptiva es una reescritura de los datos, y no hemos resumido nada todavía. Pero es apropiada para empezar a contestar la pregunta:\n\n¿Qué indican estos datos acerca de qué tratamiento es mejor? ¿Acerca del tamaño de cálculos grandes o chicos?\n\nSupongamos que otro analista decide comparar los pacientes que recibieron cada tratamiento, ignorando la variable de tamaño:\n\ncalculos |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\ny parece ser que el tratamiento \\(B\\) es mejor que el \\(A\\). Esta es una paradoja (un ejemplo de la paradoja de Simpson) . Si un médico no sabe que tipo de cálculos tiene el paciente, ¿entonces debería recetar \\(B\\)? ¿Si sabe debería recetar \\(A\\)? Esta discusión parece no tener mucho sentido.\nPodemos investigar por qué está pasando esto considerando la siguiente tabla, que solo examina cómo se asignó el tratamiento dependiendo del tipo de cálculos de cada paciente:\n\ncalculos |> group_by(tratamiento, tamaño) |> count() |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nn\n\n\n\n\nA\nchicos\n87\n\n\nA\ngrandes\n263\n\n\nB\nchicos\n270\n\n\nB\ngrandes\n80\n\n\n\n\n\n\n\nNuestra hipótesis aquí es que la decisión de qué tratamiento usar depende del tamaño de los cálculos. En este caso, hay una decisión pues A es una cirugía y B es un procedimiento menos invasivo, y se prefiere utilizar el tratamiento \\(A\\) para cálculos grandes, y \\(B\\) para cálculos chicos. Esto quiere decir que en la tabla total el tratamiento \\(A\\) está en desventaja porque se usa en casos más difíciles, pero el tratamiento \\(A\\) parece ser en general mejor. La razón es probablemente un proceso de optimización de recursos y riesgo que hacen los doctores.\n\nEn este caso, una mejor respuesta a la pregunta de qué tratamiento es mejor es la que presenta los datos desagregados.\nLa tabla desagregada de asignación del tratamiento nos informa acerca de cómo se está distribuyendo el tratamiento en los pacientes.\n\n\n\n\n\n\n\nNota\n\n\n\nLos resúmenes descriptivos acompañados de hipótesis causales acerca del proceso generador de datos, nos guía hacia descripciones interpretables de los datos.\n\n\nLas explicaciones no son tan simples y, otra vez, interviene el comportamiento de doctores, tratamientos, y distintos tipos de padecimientos.\nPodemos codificar la información causal con un diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n T \n M \n C\n edge [minlen = 3]\n T -> M\n C -> T\n C -> M\n{ rank = same; M; T }\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).\n\n\nEjemplo (cálculos renales 2)\nContrastemos el ejemplo anterior usando exactamente la misma tabla de datos, pero con el supuesto de un proceso generador diferente. En este caso, los tratamientos son para mejorar alguna enfermedad del corazón. Sabemos que parte del efecto de este tratamiento ocurre gracias a una baja en presión arterial de los pacientes, así que después de administrar el tratamiento, se toma la presión arterial de los pacientes. Ahora tenemos la tabla agregada y desagregada como sigue:\n\ncorazon <- calculos |> \n select(tratamiento, presión = tamaño, resultado) |> \n mutate(presión = ifelse(presión == \"grandes\", \"alta\", \"baja\"))\ncorazon_agregada <- corazon |> \n group_by(tratamiento, presión, resultado) |> \n count()\ncorazon_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, presión, total, prop_mejora) |> \n arrange(presión) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\npresión\ntotal\nprop_mejora\n\n\n\n\nA\nalta\n263\n0.73\n\n\nB\nalta\n80\n0.69\n\n\nA\nbaja\n87\n0.93\n\n\nB\nbaja\n270\n0.87\n\n\n\n\n\n\n\n\ncorazon |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\n¿Cuál creemos que es el mejor tratamiento en este caso? ¿Deberíamos usar la tabla agregada o la desagregada por presión?\n\nEn este caso, la tabla agregada es más apropiada (B es mejor tratamiento).\nLa razón es que presión en este caso es una consecuencia de tomar el tratamiento, y como las tablas muestran, B es más exitoso en bajar la presión de los pacientes.\nSi sólo comparamos dentro de los grupos de presión baja o de presión alta, ignoramos lo más importante del tratamiento en la probabilidad de mejorar.\n\nNuestros supuestos causales podemos mostrarlos con el siguiente diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n P\n T \n M \n edge [minlen = 3]\n T -> P\n P -> M\n T -> M\n{ rank = same; M; T}\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nNótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas." }, { "objectID": "01-introduccion.html#diagramas-causales-1", @@ -67,7 +67,7 @@ "href": "02-flujo-basico.html#paso-4-probar-el-proceso-de-estimación", "title": "2  Flujo de trabajo básico: motivación", "section": "2.4 Paso 4: Probar el proceso de estimación", - "text": "2.4 Paso 4: Probar el proceso de estimación\nAntes de utilizar datos, verificamos cómo se comporta nuestro proceso de estimación de acuerdo a los supuestos de nuestro modelo generativo.\n\n\n\n\n\n\nVerificación a priori\n\n\n\nLo 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, y que las estimaciones que arroja son apropiadas para la cantidad que nos interesa estimar. El procedimiento a grandes rasgos es:\n\nEstablecer valores de los parámetros a estimar\nSimular 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).\nCalcular posterior de las cantidades de interés\nCompara los valores de 1) con la posterior de 3)\n\n\n\nDefinir 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:\n\ntheta <- 0.2\nN <- 30\n# simular\nset.seed(9914)\ndatos_sim <- sim_pos_neg(theta = theta, N = N)\nposterior <- calcular_posterior(datos_sim$Pos)\nggplot(posterior, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(xintercept = theta, color = \"red\", linetype = \"dashed\")\n\n\n\n\nEn este caso, la estimación parece correcta. Podemo repetir el proceso con distintos valores de \\(\\theta\\):\n\nset.seed(21)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular valor inicial\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n posterior <- calcular_posterior(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nY vemos que en general nuestro método parece funcionar correctamente.\n\n\n\n\n\n\nObservaciones\n\n\n\n\nMá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.\nTí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.\nVeremos 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.\n\n\n\nEste 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:\n\nsimulacion_rep <- map_df(1:20, \n function(rep){\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 3)\n posterior <- calcular_posterior(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nNuestra 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\\).\n\n2.4.1 Introduciendo un bug\nSupongamos que tenemos un error en el cálculo de la posterior:\n\ncalcular_posterior_bug <- function(muestra, prob_priori){\n # distribución inicial o a prior\n theta <- seq(0, 1, length.out = 11)\n priori <- tibble(theta = theta, prob_priori = (1 - theta) * (1 - theta)) |> \n mutate(prob_priori = prob_priori / sum(prob_priori))\n # calcular la probabilidad posterior\n N <- length(muestra)\n Npos <- sum(muestra)\n prob_post <- tibble(theta = theta) |> \n left_join(priori, by = \"theta\") |> \n # la siguiente línea tiene un error!\n mutate(prob_posterior = theta ^ Npos * (1 - theta)^((N - Npos * prob_priori))) |> \n mutate(prob_posterior = prob_posterior / sum(prob_posterior)) \n prob_post |> select(theta, prob_posterior)\n}\n\nNuestro chequeo apriori se ve entonces:\n\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular valor inicial\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n posterior <- calcular_posterior_bug(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nDonde vemos en varios casos que la “posterior” está lejos de ser consistente con los valores simulados de prueba para \\(\\theta\\).\n\n\n\n\n\n\nAspectos numéricos\n\n\n\nEs 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.\n\n\nAunque en este caso no es crítico, la siguiente función sigue esta práctica que en general es necesario seguir:\n\n# Evitar desbordes al sumar exponenciales\nlog_sum_exp <- function(x){\n max_x <- max(x)\n max_x + log(sum(exp(x - max_x)))\n}\ncalcular_posterior <- function(muestra, prob_priori){\n # evitar 0 o 1 exactos\n theta <- seq(1e-12, 1 - 1e-12, length.out = 11)\n # no es necesario normalizar esta distribución apriori\n log_priori <- tibble(theta = theta, log_prob_priori = 2 * log(1 - theta)) \n # calcular la probabilidad posterior\n N <- length(muestra)\n Npos <- sum(muestra)\n prob_post_tbl <- tibble(theta = theta) |> \n left_join(log_priori, by = \"theta\") |> \n # log verosimilitud\n mutate(log_prob_posterior = \n Npos * log(theta) + log(1 - theta) * (N - Npos)) |> \n # sumar log apriori\n mutate(log_prob_posterior = log_prob_posterior + log_prob_priori) |> \n mutate(log_prob_posterior_norm = \n log_prob_posterior - log_sum_exp(log_prob_posterior)) |> \n mutate(prob_posterior = exp(log_prob_posterior_norm))\n prob_post_tbl |> select(theta, prob_posterior)\n}\n\nEjercicio: 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." + "text": "2.4 Paso 4: Probar el proceso de estimación\nAntes de utilizar datos, verificamos cómo se comporta nuestro proceso de estimación de acuerdo a los supuestos de nuestro modelo generativo.\n\n\n\n\n\n\nVerificación a priori\n\n\n\nLo 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, y que las estimaciones que arroja son apropiadas para la cantidad que nos interesa estimar. El procedimiento a grandes rasgos es:\n\nEstablecer valores de los parámetros a estimar\nSimular 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).\nCalcular posterior de las cantidades de interés\nCompara los valores de 1) con la posterior de 3)\n\n\n\nDefinir 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:\n\ntheta <- 0.2\nN <- 30\n# simular\nset.seed(9914)\ndatos_sim <- sim_pos_neg(theta = theta, N = N)\nposterior <- calcular_posterior(datos_sim$Pos)\nggplot(posterior, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(xintercept = theta, color = \"red\", linetype = \"dashed\")\n\n\n\n\nEn este caso, la estimación parece correcta. Podemo repetir el proceso con distintos valores de \\(\\theta\\):\n\nset.seed(21)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular valor inicial\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n posterior <- calcular_posterior(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nY vemos que en general nuestro método parece funcionar correctamente.\n\n\n\n\n\n\nObservaciones\n\n\n\n\nMá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.\nTí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.\nVeremos más de chequeos o pruebas predictivas a priori, que en general también sirven para entender la adecuación del modelo y supuestos en términos de teoría y\n\n\n\nEste 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:\n\nsimulacion_rep <- map_df(1:20, \n function(rep){\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 3)\n posterior <- calcular_posterior(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nNuestra 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\\).\n\n2.4.1 Introduciendo un bug\nSupongamos que tenemos un error en el cálculo de la posterior:\n\ncalcular_posterior_bug <- function(muestra, prob_priori){\n # distribución inicial o a prior\n theta <- seq(0, 1, length.out = 11)\n priori <- tibble(theta = theta, prob_priori = (1 - theta) * (1 - theta)) |> \n mutate(prob_priori = prob_priori / sum(prob_priori))\n # calcular la probabilidad posterior\n N <- length(muestra)\n Npos <- sum(muestra)\n prob_post <- tibble(theta = theta) |> \n left_join(priori, by = \"theta\") |> \n # la siguiente línea tiene un error!\n mutate(prob_posterior = theta ^ Npos * (1 - theta)^((N - Npos * prob_priori))) |> \n mutate(prob_posterior = prob_posterior / sum(prob_posterior)) \n prob_post |> select(theta, prob_posterior)\n}\n\nNuestro chequeo apriori se ve entonces:\n\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular valor inicial\n theta_sim <- sample(seq(0, 1, length.out = 11), \n prob = prob_priori$prob_priori, size = 1)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n posterior <- calcular_posterior_bug(datos_sim$Pos)\n posterior |> mutate(theta = theta) |> \n mutate(rep = rep) |> \n mutate(theta_sim = theta_sim)\n })\nggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nDonde vemos en varios casos que la “posterior” está lejos de ser consistente con los valores simulados de prueba para \\(\\theta\\).\n\n\n\n\n\n\nAspectos numéricos\n\n\n\nEs 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.\n\n\nAunque en este caso no es crítico, la siguiente función sigue esta práctica que en general es necesario seguir:\n\n# Evitar desbordes al sumar exponenciales\nlog_sum_exp <- function(x){\n max_x <- max(x)\n max_x + log(sum(exp(x - max_x)))\n}\ncalcular_posterior <- function(muestra, prob_priori){\n # evitar 0 o 1 exactos\n theta <- seq(1e-12, 1 - 1e-12, length.out = 11)\n # no es necesario normalizar esta distribución apriori\n log_priori <- tibble(theta = theta, log_prob_priori = 2 * log(1 - theta)) \n # calcular la probabilidad posterior\n N <- length(muestra)\n Npos <- sum(muestra)\n prob_post_tbl <- tibble(theta = theta) |> \n left_join(log_priori, by = \"theta\") |> \n # log verosimilitud\n mutate(log_prob_posterior = \n Npos * log(theta) + log(1 - theta) * (N - Npos)) |> \n # sumar log apriori\n mutate(log_prob_posterior = log_prob_posterior + log_prob_priori) |> \n mutate(log_prob_posterior_norm = \n log_prob_posterior - log_sum_exp(log_prob_posterior)) |> \n mutate(prob_posterior = exp(log_prob_posterior_norm))\n prob_post_tbl |> select(theta, prob_posterior)\n}\n\nEjercicio: 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." }, { "objectID": "02-flujo-basico.html#paso-5-analizar-los-datos-y-resumir-resultados.", @@ -88,7 +88,7 @@ "href": "02-flujo-basico.html#versión-continua", "title": "2  Flujo de trabajo básico: motivación", "section": "2.7 Versión continua", - "text": "2.7 Versión continua\nEn 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]\\).\nPara 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:\n\\[ 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:\n\\[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\n\\[p(\\theta|D) \\propto \\theta^{N_+}(1-\\theta)^{N_-}\\]\nEn este caso, para normalizar tenemos que hacer la integral de la expresión de la derecha, y dividir por el resultado. En general, escribiremos\n\\[B(a,b) = \\int_{0}^1 \\theta^{a-1}(1-\\theta)^{b-1} d\\theta\\] así que en nuestro caso, la posterior es:\n\\[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\\).\nPor ejemplo, si observamos 2 positivos y tres negativos, nuestra posterior es una beta con parámetros \\((3,4)\\), y se ve así:\n\nlibrary(tidyverse)\ntheta <- seq(0,1, 0.01)\ntibble(theta = theta, densidad = dbeta(theta, 3, 4)) |> \n ggplot(aes(x = theta, y = densidad)) +\n geom_line() +\n labs(x = \"theta\", y = \"Densidad posterior\") \n\n\n\n\nNotamos 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\n\\[p(\\theta) \\propto \\theta^{a}(1-\\theta)^{b}\\]\nentonces la posterior, por la fórmula de Bayes, es:\n\\[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)\\).\n\n2.7.1 Ejercicio: actualizaciones de posterior\nPodemos examinar la posterior para dados distintos datos. Supondremos que la distribución a priori es uniforme.\n\nset.seed(92192)\ntheta_seq <- seq(0,1, 0.001)\ndatos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> \n mutate(obs = ifelse(Pos==1, \"P\", \"N\")) |> \n mutate(n = 1:12)\n# graficar posteriores para cada n\ndatos_graf <- datos_sim |> \n mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> \n mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> \n group_by(n) |>\n mutate(dens_graf = \n list(tibble(theta = theta_seq, densidad = dbeta(theta_seq, 1 + n_pos, 1 + n_neg)))) |> \n unnest(dens_graf)\nggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +\n geom_line() + \n facet_wrap(~ muestra) +\n geom_abline(slope = 0, intercept = 1, color = \"gray\") \n\n\n\n\nAhora 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)\\):\n\nset.seed(92192)\ntheta_seq <- seq(0,1, 0.001)\ndatos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> \n mutate(obs = ifelse(Pos==1, \"P\", \"N\")) |> \n mutate(n = 1:12)\n# graficar posteriores para cada n\ndatos_graf <- datos_sim |> \n mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> \n mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> \n group_by(n) |>\n mutate(dens_graf = \n list(tibble(theta = theta_seq, \n densidad = dbeta(theta_seq, 1 + n_pos + 0, 1 + n_neg + 2)))) |> \n unnest(dens_graf)\nggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +\n geom_line() + \n facet_wrap(~ muestra) +\n geom_abline(slope = -2, intercept = 2, color = \"gray\") \n\n\n\n\n\nEn 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).\nAunque 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.\n\n\n2.7.2 Métodos Monte Carlo\nUna 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 consisten de integrales, por ejemplo:\n\nLa media o mediana posterior\nDeciles o u otro tipo de percentiles de la posterior\nIntervalos de probabilidad posterior\n\nEste 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 lugar de usar una función de calcular_posterior, construimos una que es simular_posterior.\nEsta función será simple porque simular de una beta es un problema estándar, y existen muchas implementaciones. Podríamos escribir, por ejemplo:\n\nsimular_posterior <- function(muestra, n){\n tibble(theta = \n rbeta(n, sum(muestra) + 1, length(muestra) - sum(muestra) + 1))\n}\n\n\nmuestra\n\n[1] 1 0 0 1 0\n\nsims_post <- simular_posterior(muestra, 10000)\n\n\nsims_post |> \n ggplot(aes(x = theta)) +\n geom_histogram(bins = 50)\n\n\n\n\nSi queremos calcular la media, por ejemplo, hacemos\n\n sims_post |> pull(theta) |> mean()\n\n[1] 0.4307996\n\n\nSi queremos la probabilidad de que la prevalencia esté por debajo de 20% hacemos:\n\nsims_post |> \n summarise(prob = mean(theta < 0.2))\n\n# A tibble: 1 × 1\n prob\n <dbl>\n1 0.0893\n\n\n\n\n\n\n\n\nMétodos Monte Carlo\n\n\n\nLos 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):\n\\[\\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\\).\n\n\nNota 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.\nFinalmente, 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:\n\nGeneramos un valor de la apriori \\(\\theta_{sim} \\sim \\text{Beta}(1,3)\\)\nSimulamos datos de la muestra (\\(N=25\\)) con el valor simulado de \\(\\theta\\)\nSimulamos 10000 valores de la posterior\nRepetimos 1-3\n\n\nset.seed(812)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 25)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior(datos_sim$Pos, 10000)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)" + "text": "2.7 Versión continua\nEn 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]\\).\nPara 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:\n\\[ 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:\n\\[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\n\\[p(\\theta|D) \\propto \\theta^{N_+}(1-\\theta)^{N_-}\\]\nEn este caso, para normalizar tenemos que hacer la integral de la expresión de la derecha, y dividir por el resultado. En general, escribiremos\n\\[B(a,b) = \\int_{0}^1 \\theta^{a-1}(1-\\theta)^{b-1} d\\theta\\] así que en nuestro caso, la posterior es:\n\\[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\\).\nPor ejemplo, si observamos 2 positivos y tres negativos, nuestra posterior es una beta con parámetros \\((3,4)\\), y se ve así:\n\nlibrary(tidyverse)\ntheta <- seq(0,1, 0.01)\ntibble(theta = theta, densidad = dbeta(theta, 3, 4)) |> \n ggplot(aes(x = theta, y = densidad)) +\n geom_line() +\n labs(x = \"theta\", y = \"Densidad posterior\") \n\n\n\n\nNotamos 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\n\\[p(\\theta) \\propto \\theta^{a}(1-\\theta)^{b}\\]\nentonces la posterior, por la fórmula de Bayes, es:\n\\[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)\\).\n\n2.7.1 Ejercicio: actualizaciones de posterior\nPodemos examinar la posterior para dados distintos datos. Supondremos que la distribución a priori es uniforme.\n\nset.seed(92192)\ntheta_seq <- seq(0,1, 0.001)\ndatos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> \n mutate(obs = ifelse(Pos==1, \"P\", \"N\")) |> \n mutate(n = 1:12)\n# graficar posteriores para cada n\ndatos_graf <- datos_sim |> \n mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> \n mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> \n group_by(n) |>\n mutate(dens_graf = \n list(tibble(theta = theta_seq, densidad = dbeta(theta_seq, 1 + n_pos, 1 + n_neg)))) |> \n unnest(dens_graf)\nggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +\n geom_line() + \n facet_wrap(~ muestra) +\n geom_abline(slope = 0, intercept = 1, color = \"gray\") \n\n\n\n\nAhora 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)\\):\n\nset.seed(92192)\ntheta_seq <- seq(0,1, 0.001)\ndatos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> \n mutate(obs = ifelse(Pos==1, \"P\", \"N\")) |> \n mutate(n = 1:12)\n# graficar posteriores para cada n\ndatos_graf <- datos_sim |> \n mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> \n mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> \n group_by(n) |>\n mutate(dens_graf = \n list(tibble(theta = theta_seq, \n densidad = dbeta(theta_seq, 1 + n_pos + 0, 1 + n_neg + 2)))) |> \n unnest(dens_graf)\nggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +\n geom_line() + \n facet_wrap(~ muestra) +\n geom_abline(slope = -2, intercept = 2, color = \"gray\") \n\n\n\n\n\nEn 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 de ser menor que un valor, o cómo se calcula la media).\nAunque 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.\n\n\n2.7.2 Métodos Monte Carlo\nUna 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 consisten de integrales, por ejemplo:\n\nLa media o mediana posterior\nDeciles o u otro tipo de percentiles de la posterior\nIntervalos de probabilidad posterior\n\nEste 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 lugar de usar una función de calcular_posterior, construimos una que es simular_posterior.\nEsta función será simple porque simular de una beta es un problema estándar, y existen muchas implementaciones. Podríamos escribir, por ejemplo:\n\nsimular_posterior <- function(muestra, n){\n tibble(theta = \n rbeta(n, sum(muestra) + 1, length(muestra) - sum(muestra) + 1))\n}\n\n\nmuestra\n\n[1] 1 0 0 1 0\n\nsims_post <- simular_posterior(muestra, 10000)\n\n\nsims_post |> \n ggplot(aes(x = theta)) +\n geom_histogram(bins = 50)\n\n\n\n\nSi queremos calcular la media, por ejemplo, hacemos\n\n sims_post |> pull(theta) |> mean()\n\n[1] 0.4307996\n\n\nSi queremos la probabilidad de que la prevalencia esté por debajo de 20% hacemos:\n\nsims_post |> \n summarise(prob = mean(theta < 0.2))\n\n# A tibble: 1 × 1\n prob\n <dbl>\n1 0.0893\n\n\n\n\n\n\n\n\nMétodos Monte Carlo\n\n\n\nLos 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):\n\\[\\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\\).\n\n\nNota 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.\nNota 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:\n\\[\\theta_i \\sim p(\\theta|D) \\implies \\frac{1}{M} \\sum_{i=1}^{M} f(\\theta_i) \\to I.\\]\n\nFinalmente, 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:\n\nGeneramos un valor de la apriori \\(\\theta_{sim} \\sim \\text{Beta}(1,3)\\)\nSimulamos datos de la muestra (\\(N=25\\)) con el valor simulado de \\(\\theta\\)\nSimulamos 10000 valores de la posterior\nRepetimos 1-3\n\n\nset.seed(812)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 25)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior(datos_sim$Pos, 10000)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)" }, { "objectID": "02-flujo-basico.html#observaciones-1", @@ -102,13 +102,13 @@ "href": "02-flujo-basico-2.html#prevalencia-con-error-conocido", "title": "3  Flujo de trabajo básico: refinando el modelo", "section": "3.1 Prevalencia con error conocido", - "text": "3.1 Prevalencia con error conocido\nNuestro 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:\n\nEn 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.\nSin considerar la incertidumbre, esto implica que la prueba tiene una sensibilidad de 84% y una especificidad de 99.5%.\n\n\n3.1.1 Paso 1: modelo generativo\nPrimero supondremos que estos porcentajes de error son fijos. Nuestro modelo que incluye el error de medición se como sigue:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n p\n Npos\n node [shape=plaintext]\n N\n Npos [label = <N<SUB>+</SUB>>]\n Nobs\n #Nneg [label = <N<SUB>-</SUB>>]\n #sens\n #esp\n edge [minlen = 3]\n p -> Npos\n #p -> Nneg\n N -> Npos\n Npos -> Nobs\n #N -> Nneg\n esp -> Nobs\n sens -> Nobs\n #esp -> Nneg\n #sens -> Nneg\n{ rank = same; p; N }\n{ rank = same; Npos}\n{ rank = max; sens; esp}\n}\n\")#, width = 200, height = 50)\n\n\n\n\n\n\nDonde 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.\nY 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:\n\nsim_pos_neg <- function(theta = 0.01, N = 20, sens = 0.84, esp = 0.995) {\n # verdaderos positivos que capturamos en la muestra\n Pos_verdadero <- rbinom(N, 1, theta)\n Neg_verdadero <- 1 - Pos_verdadero\n # positivos observados en la muestra: si es positivo, calculamos\n # la probabilidad de que realmente sea positivo\n sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> \n mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> \n mutate(Neg = 1 - Pos)\n # Observaciones\n sim_tbl |> select(Pos, Neg)\n}\n\nHacemos unas pruebas:\n\nset.seed(8212)\nsim_pos_neg(theta = 0.3, N = 1e7, sens = 0.7, esp = 1) |> pull(Pos) |> mean() |> \n round(4)\n\n[1] 0.2099\n\nsim_pos_neg(theta = 0.3, N = 1e7, sens = 1, esp = 1) |> pull(Pos) |> mean() |> \n round(4)\n\n[1] 0.3001\n\n\n\n\n3.1.2 Paso 2: cantidad a estimar\nEn 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.\n\n\n3.1.3 Paso 3: modelo estadístico\nEl 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:\n\\(\\theta_{obs} = P(Positivo | \\theta, sens, esp) = \\theta \\cdot sens + (1-\\theta) \\cdot (1-esp)\\)\nSi llamamos a esta cantidad \\(\\theta_{obs}\\), de forma que dada una muestra de 0’s y 1’s, tenemos\n\\[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:\n\\[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:\n\\[p(\\theta|D, sens, esp) \\propto (\\theta \\cdot sens + (1-\\theta) \\cdot (1-esp))^{N_{+}}(\\theta(1-sens) + (1-\\theta)esp)^{N_{-}}\\]\nEsta 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:\n\n\n\n\n\n\nMétodo de aproximación de rejilla\n\n\n\n\nDividimos 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\nNormalizamos estos valores para que sumen 1, y obtenemos una distribución discreta que aproxima la posterior\nMuestreamos de esta distribución discreta para obtener una muestra de la posterior\n\nEste 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.\n\n\nAquí implementamos esta técnica de aproximación por rejilla. Incluimos también una Beta(1,3) como a priori:\n\nsimular_posterior_error <- function(muestra, n, sens = 1, esp = 1){\n theta <- seq(1e-12, 1-1e-12, by = 0.0001)\n p_obs <- theta * sens + (1 - theta) * (1 - esp)\n # verosimilitud (en logaritmo)\n log_dens_sin_norm <- log(p_obs) * sum(muestra) + \n log(1-p_obs) * (length(muestra) - sum(muestra))\n # a priori\n log_dens_sin_norm <- log_dens_sin_norm + dbeta(theta, 1, 3, log = TRUE)\n # normalizar\n log_dens_norm <- log_dens_sin_norm - log_sum_exp(log_dens_sin_norm)\n densidad_post <- exp(log_dens_norm)\n tibble(theta = sample(theta, size = n, replace = TRUE, prob = densidad_post))\n}\n\nY ahora podemos ver cómo se ve la posterior:\n\nset.seed(328)\nuna_muestra <- sim_pos_neg(theta = 0.2, N = 600, sens = 0.6, esp = 0.999)\nmean(una_muestra$Pos)\n\n[1] 0.1233333\n\nsims_post_error <- \n simular_posterior_error(una_muestra$Pos, 5000, sens = 0.6, esp = 0.999) \nsims_post_error |>\n ggplot(aes(x = theta)) +\n geom_histogram()\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\nAhora seguimos el flujo. Agregaremos la verificación a priori para entender si nuestro modelo recupera los parámetros.\n\nset.seed(8112)\nsimulacion_rep_error <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior_error(datos_sim$Pos, 10000, sens = 0.6, esp = 0.999)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep_error, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nContrasta con lo que pasaría si usaramos el modelo sin considerar fuentes de error:\n\nset.seed(812)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior_error(datos_sim$Pos, 10000, 1, 1)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nEste resultado está lejos de ser aceptable.\nComparamos esta densidad con lo que obtendríamos sin considerar el error de medición, con los mismos datos:\n\nset.seed(8)\nsims_post <- \n simular_posterior_error(una_muestra$Pos, 5000, 1, 1)\nambas_sims_tbl <- \n sims_post_error |>\n mutate(tipo = \"Con error\") |>\n bind_rows(sims_post |>\n mutate(tipo = \"Sin error\"))\nambas_sims_tbl |> ggplot(aes(x = theta, fill = tipo)) +\n geom_histogram(position = \"identity\", alpha = 0.5, bins = 50) +\n scale_fill_manual(values = c(\"red\", \"blue\")) +\n geom_vline(xintercept = 0.2, linetype = \"dashed\", color = \"black\")\n\n\n\n\nY 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.\nAunque este ejemplo es claro, prevalencia, sensibilidad y especificidad interactúan de maneras a veces poco intuitivas." + "text": "3.1 Prevalencia con error conocido\nNuestro 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:\n\nEn 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.\nSin considerar la incertidumbre, esto implica que la prueba tiene una sensibilidad de 84% y una especificidad de 99.5%.\n\n\n3.1.1 Paso 1: modelo generativo\nPrimero supondremos que estos porcentajes de error son fijos. Nuestro modelo que incluye el error de medición se como sigue:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n p\n Npos\n node [shape=plaintext]\n N\n Npos [label = <N<SUB>+</SUB>>]\n Nobs [label = <N<SUB>obs</SUB>>]\n #Nneg [label = <N<SUB>-</SUB>>]\n #sens\n #esp\n edge [minlen = 3]\n p -> Npos\n #p -> Nneg\n N -> Npos\n Npos -> Nobs\n #N -> Nneg\n esp -> Nobs\n sens -> Nobs\n #esp -> Nneg\n #sens -> Nneg\n{ rank = same; p; N }\n{ rank = same; Npos}\n{ rank = max; sens; esp}\n}\n\")#, width = 200, height = 50)\n\n\n\n\n\n\nDonde 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.\nY 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:\n\nsim_pos_neg <- function(theta = 0.01, N = 20, sens = 0.84, esp = 0.995) {\n # verdaderos positivos que capturamos en la muestra\n Pos_verdadero <- rbinom(N, 1, theta)\n Neg_verdadero <- 1 - Pos_verdadero\n # positivos observados en la muestra: si es positivo, calculamos\n # la probabilidad de que realmente sea positivo\n sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> \n mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> \n mutate(Neg = 1 - Pos)\n # Observaciones\n sim_tbl |> select(Pos, Neg)\n}\n\nHacemos unas pruebas:\n\nset.seed(8212)\nsim_pos_neg(theta = 0.3, N = 1e7, sens = 0.7, esp = 1) |> pull(Pos) |> mean() |> \n round(4)\n\n[1] 0.2099\n\nsim_pos_neg(theta = 0.3, N = 1e7, sens = 1, esp = 1) |> pull(Pos) |> mean() |> \n round(4)\n\n[1] 0.3001\n\n\n\n\n3.1.2 Paso 2: cantidad a estimar\nEn 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.\n\n\n3.1.3 Paso 3: modelo estadístico\nEl 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:\n\\(\\theta_{obs} = P(Positivo | \\theta, sens, esp) = \\theta \\cdot sens + (1-\\theta) \\cdot (1-esp)\\)\nSi llamamos a esta cantidad \\(\\theta_{obs}\\), de forma que dada una muestra de 0’s y 1’s, tenemos que la verosimilitud de la muestra dada cada conjetura \\(\\theta\\), y con \\(sens\\) y \\(esp\\) fijas, es:\n\\[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:\n\\[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:\n\\[p(\\theta|D, sens, esp) \\propto (\\theta \\cdot sens + (1-\\theta) \\cdot (1-esp))^{N_{+}}(\\theta(1-sens) + (1-\\theta)esp)^{N_{-}}\\]\nEsta 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:\n\n\n\n\n\n\nMétodo de aproximación de rejilla\n\n\n\n\nDividimos 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\nNormalizamos estos valores para que sumen 1, y obtenemos una distribución discreta que aproxima la posterior\nMuestreamos de esta distribución discreta para obtener una muestra de la posterior\n\nEste 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.\n\n\nAquí implementamos esta técnica de aproximación por rejilla. Incluimos también una Beta(1,3) como a priori:\n\nsimular_posterior_error <- function(muestra, n, sens = 1, esp = 1){\n theta <- seq(1e-12, 1-1e-12, by = 0.0001)\n p_obs <- theta * sens + (1 - theta) * (1 - esp)\n # verosimilitud (en logaritmo)\n log_dens_sin_norm <- log(p_obs) * sum(muestra) + \n log(1-p_obs) * (length(muestra) - sum(muestra))\n # a priori\n log_dens_sin_norm <- log_dens_sin_norm + dbeta(theta, 1, 3, log = TRUE)\n # normalizar\n log_dens_norm <- log_dens_sin_norm - log_sum_exp(log_dens_sin_norm)\n densidad_post <- exp(log_dens_norm)\n tibble(theta = sample(theta, size = n, replace = TRUE, prob = densidad_post))\n}\n\nY ahora podemos ver cómo se ve la posterior:\n\nset.seed(328)\nuna_muestra <- sim_pos_neg(theta = 0.2, N = 600, sens = 0.6, esp = 0.999)\nmean(una_muestra$Pos)\n\n[1] 0.1233333\n\nsims_post_error <- \n simular_posterior_error(una_muestra$Pos, 5000, sens = 0.6, esp = 0.999) \nsims_post_error |>\n ggplot(aes(x = theta)) +\n geom_histogram()\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\nAhora seguimos el flujo. Agregaremos la verificación a priori para entender si nuestro modelo recupera los parámetros.\n\nset.seed(8112)\nsimulacion_rep_error <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior_error(datos_sim$Pos, 10000, sens = 0.6, esp = 0.999)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep_error, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nFigura 3.1: Verificación a priori\n\n\n\n\nContrasta con lo que pasaría si usaramos el modelo sin considerar fuentes de error:\n\nset.seed(812)\nsimulacion_rep <- map_df(1:20, \n function(rep){\n # simular de la apriori\n theta_sim <- rbeta(1, 1, 3)\n # simular datos según modelo\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999)\n # simulaciones montecarlo para la posterior\n posterior <- simular_posterior_error(datos_sim$Pos, 10000, 1, 1)\n # junta todo\n posterior |> mutate(n_sim = n()) |>\n mutate(rep = rep) |>\n mutate(theta_sim = theta_sim)\n })\n\nAhora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:\n\nggplot(simulacion_rep, aes(x = theta)) +\n geom_histogram(bins = 50) +\n labs(x = \"theta\", y = \"Prob posterior\") +\n geom_vline(aes(xintercept = theta_sim), color = \"red\", linetype = \"dashed\") +\n facet_wrap(~rep)\n\n\n\n\nFigura 3.2: Verificación a priori fallida (modelo incorrecto)\n\n\n\n\nEste resultado está lejos de ser aceptable.\nComparamos esta densidad con lo que obtendríamos sin considerar el error de medición, con los mismos datos:\n\nset.seed(8)\nsims_post <- \n simular_posterior_error(una_muestra$Pos, 5000, 1, 1)\nambas_sims_tbl <- \n sims_post_error |>\n mutate(tipo = \"Con error\") |>\n bind_rows(sims_post |>\n mutate(tipo = \"Sin error\"))\nambas_sims_tbl |> ggplot(aes(x = theta, fill = tipo)) +\n geom_histogram(position = \"identity\", alpha = 0.5, bins = 50) +\n scale_fill_manual(values = c(\"red\", \"blue\")) +\n geom_vline(xintercept = 0.2, linetype = \"dashed\", color = \"black\")\n\n\n\n\nY 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.\nAunque este ejemplo es claro, prevalencia, sensibilidad y especificidad interactúan de maneras a veces poco intuitivas." }, { "objectID": "02-flujo-basico-2.html#prevalencia-con-datos-de-referencia", "href": "02-flujo-basico-2.html#prevalencia-con-datos-de-referencia", "title": "3  Flujo de trabajo básico: refinando el modelo", "section": "3.2 Prevalencia con datos de referencia", - "text": "3.2 Prevalencia con datos de referencia\nAhora 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).\nDenotamos como \\(Ref\\) a los datos de referencia de “estándar de oro”.\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n theta\n esp\n sens\n Npos [label = <N<SUB>+</SUB>>]\n node [shape=plaintext]\n Nobs [label = <N<SUB>obs</SUB>>]\n # Nneg [label = <N<SUB>-</SUB>>]\n edge [minlen = 3]\n theta -> Npos\n #p -> Nneg\n N -> Npos\n Npos -> Nobs\n #N -> Nneg\n esp -> Nobs\n sens -> Nobs\n #esp -> Nneg\n #sens -> Nneg\n esp -> Ref\n sens -> Ref\n{ rank = same; theta; N }\n#{ rank = same; Npos; Nneg}\n{ rank = max; sens; esp}\n}\n\")#, width = 200, height = 50)\n\n\n\n\n\n\nPor 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:\n\nsim_pos_neg <- function(p = 0.01, N = 20, pos_gold = c(103,122), neg_gold = c(399,401)) {\n # Simular especificidad y sensibilidad\n sens <- rbeta(1, pos_gold[1] + 1, pos_gold[2] - pos_gold[1] + 1)\n esp <- rbeta(1, neg_gold[1] + 1, neg_gold[2] - neg_gold[1] + 1)\n # verdaderos positivos que capturamos en la muestra\n Pos_verdadero <- rbinom(N, 1, p)\n Neg_verdadero <- 1 - Pos_verdadero\n # positivos observados en la muestra: si es positivo, calculamos\n # la probabilidad de que realmente sea positivo\n sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> \n mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> \n mutate(Neg = 1 - Pos)\n # Observaciones\n sim_tbl |> select(Pos, Neg)\n}\n\nConsiderando 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. 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.\n\nlibrary(cmdstanr)\nmod_sc <- cmdstan_model(\"./src/sclara.stan\")\nprint(mod_sc)\n\ndata {\n int<lower=0> N;\n int<lower=0> n;\n int<lower=0> kit_pos;\n int<lower=0> n_kit_pos;\n int<lower=0> kit_neg;\n int<lower=0> n_kit_neg;\n}\n\nparameters {\n real<lower=0, upper=1> p; //seroprevalencia\n real<lower=0, upper=1> sens; //sensibilidad\n real<lower=0, upper=1> esp; //especificidad\n}\n\ntransformed parameters {\n real<lower=0, upper=1> prob_pos;\n\n prob_pos = p * sens + (1 - p) * (1 - esp);\n\n}\nmodel {\n // modelo de número de positivos\n n ~ binomial(N, prob_pos);\n // modelos para resultados del kit\n kit_pos ~ binomial(n_kit_pos, sens);\n kit_neg ~ binomial(n_kit_neg, esp);\n // iniciales para cantidades no medidas\n p ~ beta(1.0, 10.0);\n sens ~ beta(2.0, 1.0);\n esp ~ beta(2.0, 1.0);\n}\n\n\n\nn <- 50\nN <- 3300\ndatos_lista <- list(N = 3300, n = 50,\n kit_pos = 103, n_kit_pos = 122,\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.6 seconds.\n\nsims <- ajuste$draws(c(\"p\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"p\"))\n\n\nresumen |> select(variable, mean, q5, q95)\n\n# A tibble: 1 × 4\n variable mean q5 q95\n <chr> <dbl> <dbl> <dbl>\n1 p 0.0104 0.00243 0.0174\n\n\nY podemos graficar la posterior de la seroprevalencia:\n\nggplot(sims, aes(x = p)) + \n geom_histogram()\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\nY 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.\nPodemos 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:\n\nggplot(sims, aes(x = esp, y = p)) + geom_point() +\n xlab(\"Especificidad del kit\") + ylab(\"Prevalencia\") + geom_smooth()\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\nLa 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.\nNó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.\nY 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." + "text": "3.2 Prevalencia con datos de referencia\nAhora 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).\nDenotamos como \\(Ref\\) a los datos de referencia de “estándar de oro”.\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n theta\n esp\n sens\n Npos [label = <N<SUB>+</SUB>>]\n node [shape=plaintext]\n Nobs [label = <N<SUB>obs</SUB>>]\n # Nneg [label = <N<SUB>-</SUB>>]\n edge [minlen = 3]\n theta -> Npos\n #p -> Nneg\n N -> Npos\n Npos -> Nobs\n #N -> Nneg\n esp -> Nobs\n sens -> Nobs\n #esp -> Nneg\n #sens -> Nneg\n esp -> Ref\n sens -> Ref\n{ rank = same; theta; N }\n#{ rank = same; Npos; Nneg}\n{ rank = max; sens; esp}\n}\n\")#, width = 200, height = 50)\n\n\n\n\n\n\nPor 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:\n\nsim_pos_neg <- function(p = 0.01, N = 20, pos_gold = c(103,122), neg_gold = c(399,401)) {\n # Simular especificidad y sensibilidad\n sens <- rbeta(1, pos_gold[1] + 1, pos_gold[2] - pos_gold[1] + 1)\n esp <- rbeta(1, neg_gold[1] + 1, neg_gold[2] - neg_gold[1] + 1)\n # verdaderos positivos que capturamos en la muestra\n Pos_verdadero <- rbinom(N, 1, p)\n Neg_verdadero <- 1 - Pos_verdadero\n # positivos observados en la muestra: si es positivo, calculamos\n # la probabilidad de que realmente sea positivo\n sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> \n mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> \n mutate(Neg = 1 - Pos)\n # Observaciones\n sim_tbl |> select(Pos, Neg)\n}\n\nConsiderando 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. 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.\n\nlibrary(cmdstanr)\nmod_sc <- cmdstan_model(\"./src/sclara.stan\")\nprint(mod_sc)\n\ndata {\n int<lower=0> N;\n int<lower=0> n;\n int<lower=0> kit_pos;\n int<lower=0> n_kit_pos;\n int<lower=0> kit_neg;\n int<lower=0> n_kit_neg;\n}\n\nparameters {\n real<lower=0, upper=1> theta; //seroprevalencia\n real<lower=0, upper=1> sens; //sensibilidad\n real<lower=0, upper=1> esp; //especificidad\n}\n\ntransformed parameters {\n real<lower=0, upper=1> prob_pos;\n\n prob_pos = theta * sens + (1 - theta) * (1 - esp);\n\n}\nmodel {\n // modelo de número de positivos\n n ~ binomial(N, prob_pos);\n // modelos para resultados del kit\n kit_pos ~ binomial(n_kit_pos, sens);\n kit_neg ~ binomial(n_kit_neg, esp);\n // iniciales para cantidades no medidas\n theta ~ beta(1.0, 10.0);\n sens ~ beta(2.0, 1.0);\n esp ~ beta(2.0, 1.0);\n}\n\n\n\nn <- 50\nN <- 3300\ndatos_lista <- list(N = 3300, n = 50,\n kit_pos = 103, n_kit_pos = 122,\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.6 seconds.\n\nsims <- ajuste$draws(c(\"theta\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"theta\"))\n\n\nresumen |> select(variable, mean, q5, q95)\n\n# A tibble: 1 × 4\n variable mean q5 q95\n <chr> <dbl> <dbl> <dbl>\n1 theta 0.0104 0.00243 0.0174\n\n\nY podemos graficar la posterior de la seroprevalencia:\n\nggplot(sims, aes(x = theta)) + \n geom_histogram()\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\nY 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.\nPodemos 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:\n\nggplot(sims, aes(x = esp, y = theta)) + geom_point() +\n xlab(\"Especificidad del kit\") + ylab(\"Prevalencia\") + geom_smooth()\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\nLa 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.\nNó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.\nY notamos que aún con una muestra relativamente grande, el rango de \\(\\theta\\) es considerable: va desde valores cercanos a 0 hasta valores alrededor de 0.025-0.03." } ] \ No newline at end of file