diff --git a/.nojekyll b/.nojekyll index 5b605f1..464c5c1 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -a91c462a \ No newline at end of file +af710a48 \ No newline at end of file diff --git a/01-introduccion.html b/01-introduccion.html index c814f1e..3e745a1 100644 --- a/01-introduccion.html +++ b/01-introduccion.html @@ -287,27 +287,27 @@

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

Ejemp mejora -B +A grandes -mejora +sin_mejora B -chicos +grandes 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 09dd89d..4af978f 100644 --- a/02-flujo-basico-2.html +++ b/02-flujo-basico-2.html @@ -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.

@@ -475,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:

diff --git a/02-flujo-basico.html b/02-flujo-basico.html index edb4f0d..427ff8e 100644 --- a/02-flujo-basico.html +++ b/02-flujo-basico.html @@ -195,8 +195,10 @@

Tabla de contenidos

  • Más de verificaciones apriori
  • 2.7.2 Métodos Monte Carlo
    • -
  • 2.8 Observaciones
    • -
  • 2.9 Ventajas y desventajas de métodos bayesianos
    • +
  • 2.8 Observaciones +
    • @@ -266,8 +268,8 @@

    ", width = 300, height = 100)
    -
    - +
    +

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

    @@ -300,8 +302,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.

    @@ -935,14 +937,10 @@

    M 
    simulacion_datos_tbl |> 
   group_by(rep, theta_sim) |> 
-  summarise(Npos = sum(Pos)) |> 
+  summarise(Npos = sum(Pos), .groups = "drop") |> 
   ggplot(aes(x = Npos)) +
   geom_bar() +
   labs(x = "Número de positivos", y = "Frecuencia (muestras)")
    -
    -
    `summarise()` has grouped output by 'rep'. You can override using the `.groups`
-argument.
    -

    @@ -959,7 +957,7 @@

    M

    -

    Con este tipo de verificaciones podemos detectar las consecuencias de nuestros supuestos (incluyendo la elección de distribuciones a prior), así como otras decisiones de modelado (como la discretización).

    +

    Con este tipo de verificaciones podemos detectar las consecuencias de nuestros supuestos (incluyendo la elección de distribuciones a priori), así como otras decisiones de modelado (como la discretización).

    Conflictos con el conocimiento del área deben ser explorados para entenderlos y si es necesario corregir nuestros supuestos.

    @@ -990,39 +988,39 @@

    -
    theta <- 0.2
-N <- 30
-# simular
-set.seed(9914)
-datos_sim <- sim_pos_neg(theta = theta, N = N)
-posterior <- calcular_posterior(datos_sim$Pos)
-ggplot(posterior, aes(x = theta, y = prob_posterior)) +
-  geom_col() +
-  labs(x = "theta", y = "Prob posterior") +
-  geom_vline(xintercept = theta, color = "red", linetype = "dashed")
    +
    theta <- 0.2
+N <- 30
+# simular
+set.seed(9914)
+datos_sim <- sim_pos_neg(theta = theta, N = N)
+posterior <- calcular_posterior(datos_sim$Pos)
+ggplot(posterior, aes(x = theta, y = prob_posterior)) +
+  geom_col() +
+  labs(x = "theta", y = "Prob posterior") +
+  geom_vline(xintercept = theta, color = "red", linetype = "dashed")

    En este caso, la estimación parece correcta. Podemo repetir el proceso con distintos valores de \(\theta\):

    -
    set.seed(21)
-simulacion_rep <- map_df(1:20, 
-    function(rep){
-      # simular valor inicial
-      theta_sim <- sample(seq(0, 1, length.out = 11), 
-        prob = prob_priori$prob_priori, size = 1)
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
-      posterior <- calcular_posterior(datos_sim$Pos)
-      posterior |> mutate(theta = theta) |> 
-        mutate(rep = rep) |> 
-        mutate(theta_sim = theta_sim)
-    })
-ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
-  geom_col() +
-  labs(x = "theta", y = "Prob posterior") +
-  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
-  facet_wrap(~rep)
    +
    set.seed(21)
+simulacion_rep <- map_df(1:20, 
+    function(rep){
+      # simular valor inicial
+      theta_sim <- sample(seq(0, 1, length.out = 11), 
+        prob = prob_priori$prob_priori, size = 1)
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
+      posterior <- calcular_posterior(datos_sim$Pos)
+      posterior |> mutate(theta = theta) |> 
+        mutate(rep = rep) |> 
+        mutate(theta_sim = theta_sim)
+    })
+ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
+  geom_col() +
+  labs(x = "theta", y = "Prob posterior") +
+  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
+  facet_wrap(~rep)

    @@ -1041,27 +1039,27 @@

    -
    simulacion_rep <- map_df(1:20, 
-    function(rep){
-      theta_sim <- sample(seq(0, 1, length.out = 11), 
-        prob = prob_priori$prob_priori, size = 1)
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 3)
-      posterior <- calcular_posterior(datos_sim$Pos)
-      posterior |> mutate(theta = theta) |> 
-        mutate(rep = rep) |> 
-        mutate(theta_sim = theta_sim)
-    })
-ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
-  geom_col() +
-  labs(x = "theta", y = "Prob posterior") +
-  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
-  facet_wrap(~rep)
    +
    simulacion_rep <- map_df(1:20, 
+    function(rep){
+      theta_sim <- sample(seq(0, 1, length.out = 11), 
+        prob = prob_priori$prob_priori, size = 1)
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 3)
+      posterior <- calcular_posterior(datos_sim$Pos)
+      posterior |> mutate(theta = theta) |> 
+        mutate(rep = rep) |> 
+        mutate(theta_sim = theta_sim)
+    })
+ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
+  geom_col() +
+  labs(x = "theta", y = "Prob posterior") +
+  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
+  facet_wrap(~rep)

    @@ -1071,40 +1069,40 @@

    2.4.1 Introduciendo un bug

    Supongamos que tenemos un error en el cálculo de la posterior:

    -
    calcular_posterior_bug <- function(muestra, prob_priori){
-  # distribución inicial o a prior
-  theta <- seq(0, 1, length.out = 11)
-  priori <- tibble(theta = theta, prob_priori = (1 - theta) * (1 - theta)) |> 
-    mutate(prob_priori = prob_priori / sum(prob_priori))
-  # calcular la probabilidad posterior
-  N <- length(muestra)
-  Npos <- sum(muestra)
-  prob_post <- tibble(theta = theta) |> 
-      left_join(priori, by = "theta") |> 
-    # la siguiente línea tiene un error!
-      mutate(prob_posterior = theta ^ Npos * (1 - theta)^((N - Npos * prob_priori))) |> 
-    mutate(prob_posterior = prob_posterior / sum(prob_posterior)) 
-  prob_post |> select(theta, prob_posterior)
-}
    +
    calcular_posterior_bug <- function(muestra, prob_priori){
+  # distribución inicial o a prior
+  theta <- seq(0, 1, length.out = 11)
+  priori <- tibble(theta = theta, prob_priori = (1 - theta) * (1 - theta)) |> 
+    mutate(prob_priori = prob_priori / sum(prob_priori))
+  # calcular la probabilidad posterior
+  N <- length(muestra)
+  Npos <- sum(muestra)
+  prob_post <- tibble(theta = theta) |> 
+      left_join(priori, by = "theta") |> 
+    # la siguiente línea tiene un error!
+      mutate(prob_posterior = theta ^ Npos * (1 - theta)^((N - Npos * prob_priori))) |> 
+    mutate(prob_posterior = prob_posterior / sum(prob_posterior)) 
+  prob_post |> select(theta, prob_posterior)
+}

    Nuestro chequeo apriori se ve entonces:

    -
    simulacion_rep <- map_df(1:20, 
-    function(rep){
-      # simular valor inicial
-      theta_sim <- sample(seq(0, 1, length.out = 11), 
-        prob = prob_priori$prob_priori, size = 1)
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
-      posterior <- calcular_posterior_bug(datos_sim$Pos)
-      posterior |> mutate(theta = theta) |> 
-        mutate(rep = rep) |> 
-        mutate(theta_sim = theta_sim)
-    })
-ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
-  geom_col() +
-  labs(x = "theta", y = "Prob posterior") +
-  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
-  facet_wrap(~rep)
    +
    simulacion_rep <- map_df(1:20, 
+    function(rep){
+      # simular valor inicial
+      theta_sim <- sample(seq(0, 1, length.out = 11), 
+        prob = prob_priori$prob_priori, size = 1)
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
+      posterior <- calcular_posterior_bug(datos_sim$Pos)
+      posterior |> mutate(theta = theta) |> 
+        mutate(rep = rep) |> 
+        mutate(theta_sim = theta_sim)
+    })
+ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) +
+  geom_col() +
+  labs(x = "theta", y = "Prob posterior") +
+  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
+  facet_wrap(~rep)

    @@ -1125,31 +1123,31 @@

    <

    Aunque en este caso no es crítico, la siguiente función sigue esta práctica que en general es necesario seguir:

    -
    # Evitar desbordes al sumar exponenciales
-log_sum_exp <- function(x){
-  max_x <- max(x)
-  max_x + log(sum(exp(x - max_x)))
-}
-calcular_posterior <- function(muestra, prob_priori){
-  # evitar 0 o 1 exactos
-  theta <- seq(1e-12, 1 - 1e-12, length.out = 11)
-  # no es necesario normalizar esta distribución apriori
-  log_priori <- tibble(theta = theta, log_prob_priori = 2 * log(1 - theta)) 
-  # calcular la probabilidad posterior
-  N <- length(muestra)
-  Npos <- sum(muestra)
-  prob_post_tbl <- tibble(theta = theta) |> 
-    left_join(log_priori, by = "theta") |> 
-    # log verosimilitud
-    mutate(log_prob_posterior = 
-        Npos * log(theta) + log(1 - theta) * (N - Npos)) |> 
-    # sumar log apriori
-    mutate(log_prob_posterior = log_prob_posterior + log_prob_priori) |> 
-    mutate(log_prob_posterior_norm = 
-      log_prob_posterior - log_sum_exp(log_prob_posterior)) |> 
-    mutate(prob_posterior = exp(log_prob_posterior_norm))
-  prob_post_tbl |> select(theta, prob_posterior)
-}
    +
    # Evitar desbordes al sumar exponenciales
+log_sum_exp <- function(x){
+  max_x <- max(x)
+  max_x + log(sum(exp(x - max_x)))
+}
+calcular_posterior <- function(muestra, prob_priori){
+  # evitar 0 o 1 exactos
+  theta <- seq(1e-12, 1 - 1e-12, length.out = 11)
+  # no es necesario normalizar esta distribución apriori
+  log_priori <- tibble(theta = theta, log_prob_priori = 2 * log(1 - theta)) 
+  # calcular la probabilidad posterior
+  N <- length(muestra)
+  Npos <- sum(muestra)
+  prob_post_tbl <- tibble(theta = theta) |> 
+    left_join(log_priori, by = "theta") |> 
+    # log verosimilitud
+    mutate(log_prob_posterior = 
+        Npos * log(theta) + log(1 - theta) * (N - Npos)) |> 
+    # sumar log apriori
+    mutate(log_prob_posterior = log_prob_posterior + log_prob_priori) |> 
+    mutate(log_prob_posterior_norm = 
+      log_prob_posterior - log_sum_exp(log_prob_posterior)) |> 
+    mutate(prob_posterior = exp(log_prob_posterior_norm))
+  prob_post_tbl |> select(theta, prob_posterior)
+}

    Ejercicio: corre las pruebas para esta versión de la función como hicimos arriba. Este es un cambio correcto, y desde el punto de vista de desarrollo, si nuestra batería de pruebas es apropiado podemos hacerlo con más confianza.

    @@ -1159,36 +1157,36 @@

    \(N=20\) personas, con 17 negativos y 3 positivos. Calculamos la posterior:

    -
    # en nuestro modelo *no* importa el orden, verifica:
-datos_tbl <- tibble(Pos = c(rep(1, 3), rep(0, 17)))
-posterior <- calcular_posterior(muestra = datos_tbl$Pos)
-ggplot(posterior, aes(x = theta, y = prob_posterior)) +
-  geom_col() +
-  labs(x = "theta", y = "Prob posterior")
    +
    # en nuestro modelo *no* importa el orden, verifica:
+datos_tbl <- tibble(Pos = c(rep(1, 3), rep(0, 17)))
+posterior <- calcular_posterior(muestra = datos_tbl$Pos)
+ggplot(posterior, aes(x = theta, y = prob_posterior)) +
+  geom_col() +
+  labs(x = "theta", y = "Prob posterior")

    Y hay varias maneras de resumir esta posterior. Por ejemplo, podemos calcular (ojo: veremos más detalles de esto más adelante):

    -
    # Media
-posterior |> 
-  mutate(theta = theta) |> 
-  summarise(media = sum(theta * prob_posterior))
    +
    # Media
+posterior |> 
+  mutate(theta = theta) |> 
+  summarise(media = sum(theta * prob_posterior))
    # A tibble: 1 × 1
   media
   <dbl>
 1 0.166
    -
    # Intervalo de alta probabilidad 90%
-posterior |> 
-  mutate(theta = theta) |> 
-  arrange(desc(prob_posterior)) |> 
-  mutate(cumsum = cumsum(prob_posterior)) |> 
-  filter(cumsum <= 0.9) |> 
-  pull(theta) |> 
-  range()
    +
    # Intervalo de alta probabilidad 90%
+posterior |> 
+  mutate(theta = theta) |> 
+  arrange(desc(prob_posterior)) |> 
+  mutate(cumsum = cumsum(prob_posterior)) |> 
+  filter(cumsum <= 0.9) |> 
+  pull(theta) |> 
+  range()
    [1] 0.1 0.2
    @@ -1207,27 +1205,51 @@

    \[p(\theta|D) \propto \theta^{N_+}(1-\theta)^{N_-}\]

    En este caso, para normalizar tenemos que hacer la integral de la expresión de la derecha, y dividir por el resultado. En general, escribiremos

    \[B(a,b) = \int_{0}^1 \theta^{a-1}(1-\theta)^{b-1} d\theta\] así que en nuestro caso, la posterior es:

    -

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

    +

    \[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 ahora. Este tipo de densidades pertenecen a la familia beta con parámetros \((a,b)\), donde \(a>0, b>0\).

    Por ejemplo, si observamos 2 positivos y tres negativos, nuestra posterior es una beta con parámetros \((3,4)\), y se ve así:

    -
    library(tidyverse)
-theta <- seq(0,1, 0.01)
-tibble(theta = theta, densidad = dbeta(theta, 3, 4)) |> 
-  ggplot(aes(x = theta, y = densidad)) +
-  geom_line() +
-  labs(x = "theta", y = "Densidad posterior")
    +
    library(tidyverse)
+theta <- seq(0,1, 0.01)
+tibble(theta = theta, densidad = dbeta(theta, 3, 4)) |> 
+  ggplot(aes(x = theta, y = densidad)) +
+  geom_line() +
+  labs(x = "theta", y = "Densidad posterior")

    Notamos adicionalmente que es posible seleccionar otra distribución inicial que no sea la uniforme. En este caso particular es conveniente (aunque no siempre tiene sentido) usar una distribución beta, de manera que es fácil ver que si ponemos por ejemplo

    -

    \[p(\theta) \propto \theta^{a}(1-\theta)^{b}\]

    +

    \[p(\theta) \propto \theta^{a-1}(1-\theta)^{b-1}\]

    entonces la posterior, por la fórmula de Bayes, es:

    -

    \[p(\theta|D) \propto \theta^{N_+ +a }(1-\theta)^{N_{-}+b}\] que también es de la familia beta, pero con parámetros \((N_+ +a+1, N_- +b+1)\).

    +

    \[p(\theta|D) \propto \theta^{N_+ +a -1 }(1-\theta)^{N_{-}+b-1}\] que también es de la familia beta, pero con parámetros \((N_{+} +a, N_{-} +b)\).

    2.7.1 Ejercicio: actualizaciones de posterior

    Podemos examinar la posterior para dados distintos datos. Supondremos que la distribución a priori es uniforme.

    +
    set.seed(92192)
+theta_seq <- seq(0,1, 0.001)
+datos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> 
+  mutate(obs = ifelse(Pos==1, "P", "N")) |> 
+  mutate(n = 1:12)
+# graficar posteriores para cada n
+datos_graf <- datos_sim |> 
+  mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> 
+  mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> 
+  group_by(n) |>
+  mutate(dens_graf = 
+    list(tibble(theta = theta_seq, 
+      densidad = dbeta(theta_seq, n_pos + 1, n_neg + 1)))) |> 
+  unnest(dens_graf)
+ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +
+  geom_line() + 
+  facet_wrap(~ muestra) +
+  geom_abline(slope = 0, intercept = 1, color = "gray")
    +
    +

    +
    +
    +

    Ahora repetimos con una inicial beta \((0,2)\) (que equivale a observar dos negativos y ningún positivo en una muestra de 3 personas), de modo que \(p(\theta) = 2(1-\theta)\):

    +
    set.seed(92192)
 theta_seq <- seq(0,1, 0.001)
 datos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> 
@@ -1239,36 +1261,13 @@ 
      mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> 
   group_by(n) |>
   mutate(dens_graf = 
-    list(tibble(theta = theta_seq, densidad = dbeta(theta_seq, 1 + n_pos, 1 + n_neg)))) |> 
-  unnest(dens_graf)
-ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +
-  geom_line() + 
-  facet_wrap(~ muestra) +
-  geom_abline(slope = 0, intercept = 1, color = "gray")
    -
    -

    -
    -
    -

    Ahora repetimos con una inicial beta \((0,2)\) (que equivale a observar dos negativos y ningún positivo en una muestra de 3 personas), de modo que \(p(\theta) = 2(1-\theta)\):

    -
    -
    set.seed(92192)
-theta_seq <- seq(0,1, 0.001)
-datos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> 
-  mutate(obs = ifelse(Pos==1, "P", "N")) |> 
-  mutate(n = 1:12)
-# graficar posteriores para cada n
-datos_graf <- datos_sim |> 
-  mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> 
-  mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> 
-  group_by(n) |>
-  mutate(dens_graf = 
-    list(tibble(theta = theta_seq, 
-                densidad = dbeta(theta_seq, 1 + n_pos + 0, 1 + n_neg + 2)))) |> 
-  unnest(dens_graf)
-ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) +
-  geom_line() + 
-  facet_wrap(~ muestra) +
-  geom_abline(slope = -2, intercept = 2, color = "gray")
    + list(tibble(theta = theta_seq, + densidad = dbeta(theta_seq, n_pos + 1, n_neg + 3)))) |> + unnest(dens_graf) +ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) + + geom_line() + + facet_wrap(~ muestra) + + geom_abline(slope = -2, intercept = 2, color = "gray")

    @@ -1281,23 +1280,23 @@

    Más de verificaciones apriori

    Antes de continuar, sin embargo, veremos cómo se veo el chequeo predictivo a priori que consideramos en la sección de arriba.

    -
    set.seed(231)
-simulacion_datos_tbl <- map_df(1:500, 
-    function(rep){
-      # apriori seleccionada
-      theta_sim <- rbeta(1, 1, 3)
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
-      tibble(rep = rep, theta_sim = theta_sim, datos_sim)
-    })
    +
    set.seed(231)
+simulacion_datos_tbl <- map_df(1:500, 
+    function(rep){
+      # apriori seleccionada
+      theta_sim <- rbeta(1, 1, 3)
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
+      tibble(rep = rep, theta_sim = theta_sim, datos_sim)
+    })

    Podemos ver por ejemplo dónde esperamos ver el número de positivos a lo largo de distintas muestras, cuando \(N=30\):

    -
    simulacion_datos_tbl |> 
-  group_by(rep, theta_sim) |> 
-  summarise(Npos = sum(Pos)) |> 
-  ggplot(aes(x = Npos)) +
-  geom_bar() +
-  labs(x = "Número de positivos", y = "Frecuencia (muestras)")
    +
    simulacion_datos_tbl |> 
+  group_by(rep, theta_sim) |> 
+  summarise(Npos = sum(Pos)) |> 
+  ggplot(aes(x = Npos)) +
+  geom_bar() +
+  labs(x = "Número de positivos", y = "Frecuencia (muestras)") 
    `summarise()` has grouped output by 'rep'. You can override using the `.groups`
 argument.
    @@ -1309,20 +1308,20 @@

    Este resultado es consecuencia de nuestros supuestos, antes de ver los datos, y resume que esperamos con mayor probabilidad un número bajo de positivos (en una muestra de N=30), y que es muy poco probable que observemos prevalencias muy altas. Dependiendo de la situación, este puede ser un resultado aceptable.

    Un resultado no aceptable para una enfermedad que sabemos que es relativamente rara (aunque tenemos incertidumbre), por ejemplo, podría ser el siguiente:

    -
    set.seed(231)
-simulacion_datos_tbl <- map_df(1:500, 
-    function(rep){
-      # apriori seleccionada
-      theta_sim <- rbeta(1, 30, 3)
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
-      tibble(rep = rep, theta_sim = theta_sim, datos_sim)
-    })
-simulacion_datos_tbl |> 
-  group_by(rep, theta_sim) |> 
-  summarise(Npos = sum(Pos)) |> 
-  ggplot(aes(x = Npos)) +
-  geom_bar() +
-  labs(x = "Número de positivos", y = "Frecuencia (muestras)")
    +
    set.seed(231)
+simulacion_datos_tbl <- map_df(1:500, 
+    function(rep){
+      # apriori seleccionada
+      theta_sim <- rbeta(1, 30, 3)
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)
+      tibble(rep = rep, theta_sim = theta_sim, datos_sim)
+    })
+simulacion_datos_tbl |> 
+  group_by(rep, theta_sim) |> 
+  summarise(Npos = sum(Pos)) |> 
+  ggplot(aes(x = Npos)) +
+  geom_bar() +
+  labs(x = "Número de positivos", y = "Frecuencia (muestras)") 
    `summarise()` has grouped output by 'rep'. You can override using the `.groups`
 argument.
    @@ -1345,37 +1344,37 @@

    <

    Este proceso puede ser no trivial incluso para densidades posteriores conocidas. La alternativa a integrar es simular de la posterior y calcular las cantidades de interés a partir de las simulaciones. En general, esto es más fácil que integrar. En nuestro ejemplo, en lugar de usar una función de calcular_posterior, construimos una que es simular_posterior.

    Esta función será simple porque simular de una beta es un problema estándar, y existen muchas implementaciones. Podríamos escribir, por ejemplo:

    -
    simular_posterior <- function(muestra, n){
-  tibble(theta = 
-    rbeta(n, sum(muestra) + 1, length(muestra) - sum(muestra) + 1))
-}
    +
    simular_posterior <- function(muestra, n){
+  tibble(theta = 
+    rbeta(n, sum(muestra) + 1, length(muestra) - sum(muestra) + 1))
+}
    -
    muestra
    +
    muestra
    [1] 1 0 0 1 0
    -
    sims_post <- simular_posterior(muestra, 10000)
    +
    sims_post <- simular_posterior(muestra, 10000)
    -
    sims_post |> 
-  ggplot(aes(x = theta)) +
-  geom_histogram(bins = 50)
    +
    sims_post |> 
+  ggplot(aes(x = theta)) +
+  geom_histogram(bins = 50)

    Si queremos calcular la media, por ejemplo, hacemos

    -
     sims_post |> pull(theta) |>  mean()
    +
     sims_post |> pull(theta) |>  mean()
    [1] 0.4280916

    Si queremos la probabilidad de que la prevalencia esté por debajo de 20% hacemos:

    -
    sims_post |> 
-  summarise(prob = mean(theta < 0.2))
    +
    sims_post |> 
+  summarise(prob = mean(theta < 0.2))
    # A tibble: 1 × 1
     prob
@@ -1385,11 +1384,11 @@ 
    <

    Muchas veces se presentan intervalos de probabilidad posterior, por ejemplo, podríamos reportar que con 90% de probabilidad la prevalencia está en el siguiente intervalo:

    -
    sims_post |> 
-  summarise(inf = quantile(theta, 0.05),
-            sup = quantile(theta, 0.95)) |> 
-  mutate(inf = round(inf, 2),
-         sup = round(sup, 2))
    +
    sims_post |> 
+  summarise(inf = quantile(theta, 0.05),
+            sup = quantile(theta, 0.95)) |> 
+  mutate(inf = round(inf, 2),
+         sup = round(sup, 2))
    # A tibble: 1 × 2
     inf   sup
@@ -1397,7 +1396,7 @@ 
    <
 1  0.16  0.73
    -

    Observación: No hay un intervalo mágico que debe reportarse (por ejemplo 95% de probabilida es una costumbre o superstición). Hay varias maneras de construir intervalos de probabilidad. Dejaremos esta discusión para más adelante.

    +

    Observación: No hay un intervalo mágico que debe reportarse (por ejemplo 95% de probabilidad es una costumbre o superstición). Hay varias maneras de construir intervalos de probabilidad. Dejaremos esta discusión para más adelante.

    @@ -1413,39 +1412,44 @@

    <

    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:

    +

    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\) en general 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\), usaríamos \(f(\theta) = \theta\). 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. Simulamos datos de la muestra (\(N=25\)) con el valor simulado de \(\theta\)
      3. -
    3. Simulamos 10000 valores de la posterior
      4. +
    4. Simulamos un número grande \(M\) de valores de la posterior (aquí usaremos \(M=10000\))
    5. Repetimos 1-3
    -
    set.seed(812)
-simulacion_rep <- map_df(1:20, 
-    function(rep){
-      # simular de la apriori
-      theta_sim <- rbeta(1, 1, 3)
-      # simular datos según modelo
-      datos_sim <- sim_pos_neg(theta = theta_sim, N = 25)
-      # simulaciones montecarlo para la posterior
-      posterior <- simular_posterior(datos_sim$Pos, 10000)
-      # junta todo
-      posterior |> mutate(n_sim = n()) |>
-        mutate(rep = rep) |>
-        mutate(theta_sim = theta_sim)
-    })
    +
    set.seed(812)
+simulacion_rep <- map_df(1:20, 
+    function(rep){
+      # simular de la apriori
+      theta_sim <- rbeta(1, 1, 3)
+      # simular datos según modelo
+      datos_sim <- sim_pos_neg(theta = theta_sim, N = 25)
+      # simulaciones montecarlo para la posterior
+      posterior <- simular_posterior(datos_sim$Pos, 10000)
+      # junta todo
+      posterior |> mutate(n_sim = n()) |>
+        mutate(rep = rep) |>
+        mutate(theta_sim = theta_sim)
+    })
+simular_posterior <- function(muestra, n){
+  tibble(theta = 
+    rbeta(n, sum(muestra) + 1, 
+             length(muestra) - sum(muestra) + 3))
+}

    Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación:

    -
    ggplot(simulacion_rep, aes(x = theta)) +
-  geom_histogram(bins = 50) +
-  labs(x = "theta", y = "Prob posterior") +
-  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
-  facet_wrap(~rep)
    +
    ggplot(simulacion_rep, aes(x = theta)) +
+  geom_histogram(bins = 50) +
+  labs(x = "theta", y = "Prob posterior") +
+  geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") +
+  facet_wrap(~rep)

    @@ -1461,11 +1465,18 @@

    Refinaremos el proceso de checar que el cómputo (checar Monte Carlo) y la inferencia (verificación apriori) es correcta bajo nuestros supuestos.
  • Finalmente, veremos qué hacer después de hacer la estimación y que los puntos de arriba están resueltos, para tener confianza en nuestras conclusiones.
    • -

    -
    -

    2.9 Ventajas y desventajas de métodos bayesianos

    +
    +

    2.8.1 Resumen

    +

    Aquí juntamos algunas observaciones que se derivan de lo que hemos visto (flujo de trabajo y estimación bayesiana):

    +
      +

    1. Todo nuestro trabajo está fundamentado en entender qué es lo que queremos estimar dentro de un modelo generativo. Los diagramas causales nos ayudan a conectar el problema de interés con nuestros modelos, a construir modelos generativos y hacer explícitos nuestros supuestos.

      2. +

    2. El proceso de estimación siempre es el mismo: nuestro estimador es la distribución posterior, que se construye a partir de la verosimilitud y la apriori (modelo generativo). Nuestro estimador es la posterior de las cantidades de interés, que pueden resumirse de distintas maneras. Cualquier cálculo derivado de otras cantidades de interés debe considerar toda la posterior (no solo la media o la moda, etc. posterior).

      3. +

    3. Nuestro proceso incluye los chequeos predictivos a priori (basados en simulación de datos). Esto son cruciales para detectar problemas en nuestros supuestos (vs teoría) y que nuestro proceso sea internamente consistente. Esto también es una verificación de la información a priori.

      4. +

    4. Generalmente es más conveniente y práctico hacer simulaciones que calcular analíticamente la posterior o sus integrales.

      5. +
    +
    diff --git a/search.json b/search.json index 37d3a41..6ea2828 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\nA\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nA\ngrandes\nsin_mejora\n\n\nA\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nB\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\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." + "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\nchicos\nmejora\n\n\nB\ngrandes\nmejora\n\n\nA\nchicos\nmejora\n\n\nB\ngrandes\nsin_mejora\n\n\nA\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nB\ngrandes\nmejora\n\n\nB\nchicos\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", @@ -60,14 +60,14 @@ "href": "02-flujo-basico.html#paso-3-definir-un-proceso-estadístico", "title": "2  Flujo de trabajo básico: motivación", "section": "2.3 Paso 3: definir un proceso estadístico", - "text": "2.3 Paso 3: definir un proceso estadístico\nDada la información limitada que tenemos acerca de la población, esperamos tener cierta incertidumbre en nuestra estimación del valor de \\(\\theta\\). En estadística bayesiana esta incertidumbre la expresamos mediante una distribución de probabilidades sobre posibles valores del \\(\\theta\\). Si denotamos por \\(D\\) a los datos observados, nuestro objetivo es calcular o aproximar\n\\[p(\\theta|D)\\] que es una distribución sobre los posibles valores de \\(\\theta\\), una vez que tenemos información de la muestra, y que pone más masa de probabilidad sobre las conjeturas de \\(\\theta\\) que son más probables o creíbles. A esta distribución le llamamos la distribución posterior de \\(\\theta\\).\nCon esta posterior podemos hacer afirmaciones probabilísticas de la forma:\n\n¿Cuál es la probabilidad de que \\(\\theta\\) sea menor a 1%? (Muy pocos seropositivos)\n¿Cuál es la probabildad de que \\(\\theta\\) sea mayor a 80%? (Población cerca de saturación)\n\nEstas cantidades se calculan, al menos teóricamente, integrando \\(p(\\theta|D)\\) sobre los valores de \\(\\theta\\) que nos interesan, por ejemplo,\n\\[P(\\theta <= 0.01) = \\int_0^{0.01} p(\\theta|D) d\\theta\\] Nota: la integral la interpretamos como suma en el caso discreto.\nSupongamos entonces una \\(\\theta\\) dada, y que observamos la muestra \\(1,0,0,1,0\\). La probabilidad de observar esta muestra es (suponiendo observaciones independientes):\n\\[\\theta(1-\\theta)(1-\\theta)\\theta(1-\\theta) = \\theta^2(1-\\theta)^3\\] Para algunos valores de \\(\\theta\\) (posibles conjeturas acerca del valor de \\(\\theta\\)) podemos escribir una tabla como sigue (Nota: discretizamos por el momento a un número finito de valores de \\(\\theta\\) para hacer el argumento más simple):\n\ntheta <- seq(0, 1, length.out = 11)\ntibble(conjetura_theta = theta, verosimiltud = theta^2 * (1 - theta)^3) |> \n kbl(col.names = c(\"Conjetura θ\", \"p(D|θ)\"),\n escape = FALSE) \n\n\n\n\nConjetura θ\np(D|θ)\n\n\n\n\n0.0\n0.00000\n\n\n0.1\n0.00729\n\n\n0.2\n0.02048\n\n\n0.3\n0.03087\n\n\n0.4\n0.03456\n\n\n0.5\n0.03125\n\n\n0.6\n0.02304\n\n\n0.7\n0.01323\n\n\n0.8\n0.00512\n\n\n0.9\n0.00081\n\n\n1.0\n0.00000\n\n\n\n\n\n\n\nEn la tabla vemos que hay algunas conjeturas, o posibles valores de \\(\\theta\\), que tienen probabilidad considerablemente más alta que otra. La notación\n\\[p(D|\\theta)\\] significa: la probabilidad de los datos \\(D\\) dado el valor de \\(\\theta\\). Nótese que esta distribución no es la posterior que describimos arriba, y no es una distribución de probabilidad sobre \\(\\theta\\) (las probabilidades no suman uno). Esta función se llama usualmente verosimilitud de los datos, e incorpora supuestos concretos del proceso generador de los datos.\nUsando reglas de probabilidad (en particular la regla de Bayes), observamos que\n\\[p(\\theta | D) = \\frac{p(D|\\theta)p(\\theta)} { p(D)}.\\] Como \\(p(\\theta|D)\\) debe dar una distribución de probabilidad (suma o integra a 1), entonces \\(p(D)\\) debe ser una constante de normalización para el numerador de la derecha, es decir, basta escribir\n\\[p(\\theta | D) \\propto p(D|\\theta)p(\\theta) \\] Ahora es donde encontramos que tenemos que tener \\(p(\\theta)\\) para poder calcular la cantidad que nos interesa, que es la distribución posterior \\(p(\\theta|D)\\). \\(p(\\theta)\\), la distribución a priori o distribución inicial es simplemente una afirmación de dónde puede estar \\(\\theta\\), antes de observar ningún dato.\nPor el momento, podríamos poner \\(p(\\theta)\\) constante, de manera que es parte de la constante de normalización, y sólo tendríamos que normalizar como sigue:\n\ntheta <- seq(0, 1, length.out = 11)\nprob_post <- tibble(conjetura = theta, probablidad = theta^2 * (1 - theta)^3) |> \n mutate(prob_posterior = probablidad / sum(probablidad)) \nprob_post |> \n kable(col.names = c(\"Conjetura θ\", \"p(D|θ)\",\"p(θ|D)\")) |>\n kable_paper()\n\n\n\n\nConjetura θ\np(D|θ)\np(θ|D)\n\n\n\n\n0.0\n0.00000\n0.0000000\n\n\n0.1\n0.00729\n0.0437444\n\n\n0.2\n0.02048\n0.1228923\n\n\n0.3\n0.03087\n0.1852385\n\n\n0.4\n0.03456\n0.2073807\n\n\n0.5\n0.03125\n0.1875188\n\n\n0.6\n0.02304\n0.1382538\n\n\n0.7\n0.01323\n0.0793879\n\n\n0.8\n0.00512\n0.0307231\n\n\n0.9\n0.00081\n0.0048605\n\n\n1.0\n0.00000\n0.0000000\n\n\n\n\n\n\n\nCon esto, expresamos nuestro conocimiento acerca de \\(\\theta\\), después de observar los datos, con una distribución posterior de probabilidad sobre las posibles conjecturas. Este es el resultado principal de inferencia bayesiana, y es la base para tomar decisiones relativas a \\(\\theta\\).\n\nUsando información adicional\nSupongamos que tenemos información adicional acerca de \\(\\theta\\), por ejemplo, que en un experimento similar anterior alguien tomó una muestra de dos personas, y encontraron dos negativos. Tenemos entonces como creencias inciales:\n\ntheta <- seq(0, 1, length.out = 11)\nprob_priori <- tibble(conjetura = theta) |> \n mutate(prob_priori = (1 - theta) * (1 - theta)) |> \n mutate(prob_priori = prob_priori / sum(prob_priori)) \nprob_priori |>\n kable(col.names = c(\"Conjetura θ\", \"p(θ)\")) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\n\n\n\n\n0.0\n0.2597403\n\n\n0.1\n0.2103896\n\n\n0.2\n0.1662338\n\n\n0.3\n0.1272727\n\n\n0.4\n0.0935065\n\n\n0.5\n0.0649351\n\n\n0.6\n0.0415584\n\n\n0.7\n0.0233766\n\n\n0.8\n0.0103896\n\n\n0.9\n0.0025974\n\n\n1.0\n0.0000000\n\n\n\n\n\n\n\nPor ejemplo, al probabilidad inicial de que \\(\\theta\\) sea muy grande es cercana a cero, pues observamos dos negativos y ningún positivo. Ahora regresamos a considerar nuestra fórmula\n\\[p(\\theta | D) \\propto p(D|\\theta)p(\\theta), \\]\nEn este caso, la apriori o inicial tiene un efecto sobre la posterior. Reconsideramos entonces la posterior de nuestra muestra de 5 personas, y calculamos el producto de \\(P(D|\\theta)\\) por \\(p(\\theta)\\):\n\nprob_post <- prob_priori |> \n mutate(verosimilitud = theta^2 * (1 - theta)^3) |> \n mutate(prod = verosimilitud * prob_priori)\n\nprob_post|> \n kable(col.names = c(\"Conjetura θ\", \"p(θ)\", \"p(D|θ)\",\n \"p(D|θ)p(θ)\")) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\np(D|θ)\np(D|θ)p(θ)\n\n\n\n\n0.0\n0.2597403\n0.00000\n0.0000000\n\n\n0.1\n0.2103896\n0.00729\n0.0015337\n\n\n0.2\n0.1662338\n0.02048\n0.0034045\n\n\n0.3\n0.1272727\n0.03087\n0.0039289\n\n\n0.4\n0.0935065\n0.03456\n0.0032316\n\n\n0.5\n0.0649351\n0.03125\n0.0020292\n\n\n0.6\n0.0415584\n0.02304\n0.0009575\n\n\n0.7\n0.0233766\n0.01323\n0.0003093\n\n\n0.8\n0.0103896\n0.00512\n0.0000532\n\n\n0.9\n0.0025974\n0.00081\n0.0000021\n\n\n1.0\n0.0000000\n0.00000\n0.0000000\n\n\n\n\n\n\n\nY finalmente, normalizamos para encontrar la probabilidad posterior:\n\nprob_post <- prob_post |> \n mutate(prob_posterior = prod / sum(prod))\n\nprob_post|> \n kable(col.names = c(\"Conjetura θ\", \"p(θ)\", \"p(D|θ)\",\n \"p(D|θ)p(θ)\", \"p(θ|D)\"), escape = FALSE) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\np(D|θ)\np(D|θ)p(θ)\np(θ|D)\n\n\n\n\n0.0\n0.2597403\n0.00000\n0.0000000\n0.0000000\n\n\n0.1\n0.2103896\n0.00729\n0.0015337\n0.0992712\n\n\n0.2\n0.1662338\n0.02048\n0.0034045\n0.2203539\n\n\n0.3\n0.1272727\n0.03087\n0.0039289\n0.2542983\n\n\n0.4\n0.0935065\n0.03456\n0.0032316\n0.2091640\n\n\n0.5\n0.0649351\n0.03125\n0.0020292\n0.1313412\n\n\n0.6\n0.0415584\n0.02304\n0.0009575\n0.0619745\n\n\n0.7\n0.0233766\n0.01323\n0.0003093\n0.0200177\n\n\n0.8\n0.0103896\n0.00512\n0.0000532\n0.0034430\n\n\n0.9\n0.0025974\n0.00081\n0.0000021\n0.0001362\n\n\n1.0\n0.0000000\n0.00000\n0.0000000\n0.0000000\n\n\n\n\n\n\n\nLa última columna nos da el resultado final de la inferencia bayesiana. Podemos resumir algunas de sus características, por ejemplo:\n\nEs muy poco probable que la seropositividad sea mayor o igual a 0.7\nUn intervalo de 90% de probabilidad para la seropositividad es \\([0.1, 0.5]\\)\n\nLa gráfica de la posterior es:\n\nprob_post |>\n ggplot(aes(x = conjetura, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") \n\n\n\n\nAhora podemos definir, para nuestro ejemplo discretizado, la función que calcula la posterior dados los pasos 1 y 2:\n\ncalcular_posterior <- 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 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\n\nmuestra <- c(1,0,0,1,0)\n\n\ncalcular_posterior(muestra, prob_priori) \n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0.0993 \n 3 0.2 0.220 \n 4 0.3 0.254 \n 5 0.4 0.209 \n 6 0.5 0.131 \n 7 0.6 0.0620 \n 8 0.7 0.0200 \n 9 0.8 0.00344 \n10 0.9 0.000136\n11 1 0 \n\n\nProcedemos ahora a hacer algunas pruebas simples de nuestra función:\n\ncalcular_posterior(rep(0, 50)) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0.996\n 2 0.1 0.004\n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0 \n 6 0.5 0 \n 7 0.6 0 \n 8 0.7 0 \n 9 0.8 0 \n10 0.9 0 \n11 1 0 \n\ncalcular_posterior(rep(1, 50)) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0 \n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0 \n 6 0.5 0 \n 7 0.6 0 \n 8 0.7 0 \n 9 0.8 0.011\n10 0.9 0.989\n11 1 0 \n\ncalcular_posterior(c(rep(0, 100), rep(1, 100))) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0 \n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0.023\n 6 0.5 0.966\n 7 0.6 0.01 \n 8 0.7 0 \n 9 0.8 0 \n10 0.9 0 \n11 1 0 \n\n\n\n\nMás verificaciones a priori\nOtra verificación útil que podemos hacer es, una vez que hemos definido nuestro modelo generativo y un modelos estadístico asociado, generar bajo simulación datos que podríamos observar. Esto tiene como fin verificar que nuestro modelo generativo y nuestro modelo estadístico producen datos que están de acuerdo con el conocimiento experto (teoría científica o conocimiento de negocio).\nAsí que simulamos datos del modelo:\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \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 tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\n\nPodemos ver por ejemplo dónde esperamos ver el número de positivos a lo largo de distintas muestras, cuando \\(N=30\\):\n\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos)) |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n`summarise()` has grouped output by 'rep'. You can override using the `.groups`\nargument.\n\n\n\n\n\nObservamos que con nuestros supuestos, hay una probabilidad alta de observar 0 positivos (alrededor de 0.30). Esto se debe en parte a la discretización que hicimos, y que nuestra apriori pone peso considerable en prevalencia igual a cero, lo que quizá no es muy realista, y probablemente deberíamos escoger al menos una discretización más fina.\nTambién, si consideramos los supuestos como correctos, esto puede indicar el riesgo de usar una muestra chica para estimar prevalencia si esta es muy baja: es probable que obtengamos 0 observaciones positivas.\n\n\n\n\n\n\nVerificación predictiva a priori\n\n\n\nCon este tipo de verificaciones podemos detectar las consecuencias de nuestros supuestos (incluyendo la elección de distribuciones a prior), así como otras decisiones de modelado (como la discretización).\nConflictos con el conocimiento del área deben ser explorados para entenderlos y si es necesario corregir nuestros supuestos.\n\n\nEste tipo de verificaciones es muy flexible, y debe adaptarse a los aspectos del conocimiento del área que son importantes para los expertos. Podemos usar todos nuestros recursos analíticos (tablas, resúmenes, gráficas) para producir estos chequeos." + "text": "2.3 Paso 3: definir un proceso estadístico\nDada la información limitada que tenemos acerca de la población, esperamos tener cierta incertidumbre en nuestra estimación del valor de \\(\\theta\\). En estadística bayesiana esta incertidumbre la expresamos mediante una distribución de probabilidades sobre posibles valores del \\(\\theta\\). Si denotamos por \\(D\\) a los datos observados, nuestro objetivo es calcular o aproximar\n\\[p(\\theta|D)\\] que es una distribución sobre los posibles valores de \\(\\theta\\), una vez que tenemos información de la muestra, y que pone más masa de probabilidad sobre las conjeturas de \\(\\theta\\) que son más probables o creíbles. A esta distribución le llamamos la distribución posterior de \\(\\theta\\).\nCon esta posterior podemos hacer afirmaciones probabilísticas de la forma:\n\n¿Cuál es la probabilidad de que \\(\\theta\\) sea menor a 1%? (Muy pocos seropositivos)\n¿Cuál es la probabildad de que \\(\\theta\\) sea mayor a 80%? (Población cerca de saturación)\n\nEstas cantidades se calculan, al menos teóricamente, integrando \\(p(\\theta|D)\\) sobre los valores de \\(\\theta\\) que nos interesan, por ejemplo,\n\\[P(\\theta <= 0.01) = \\int_0^{0.01} p(\\theta|D) d\\theta\\] Nota: la integral la interpretamos como suma en el caso discreto.\nSupongamos entonces una \\(\\theta\\) dada, y que observamos la muestra \\(1,0,0,1,0\\). La probabilidad de observar esta muestra es (suponiendo observaciones independientes):\n\\[\\theta(1-\\theta)(1-\\theta)\\theta(1-\\theta) = \\theta^2(1-\\theta)^3\\] Para algunos valores de \\(\\theta\\) (posibles conjeturas acerca del valor de \\(\\theta\\)) podemos escribir una tabla como sigue (Nota: discretizamos por el momento a un número finito de valores de \\(\\theta\\) para hacer el argumento más simple):\n\ntheta <- seq(0, 1, length.out = 11)\ntibble(conjetura_theta = theta, verosimiltud = theta^2 * (1 - theta)^3) |> \n kbl(col.names = c(\"Conjetura θ\", \"p(D|θ)\"),\n escape = FALSE) \n\n\n\n\nConjetura θ\np(D|θ)\n\n\n\n\n0.0\n0.00000\n\n\n0.1\n0.00729\n\n\n0.2\n0.02048\n\n\n0.3\n0.03087\n\n\n0.4\n0.03456\n\n\n0.5\n0.03125\n\n\n0.6\n0.02304\n\n\n0.7\n0.01323\n\n\n0.8\n0.00512\n\n\n0.9\n0.00081\n\n\n1.0\n0.00000\n\n\n\n\n\n\n\nEn la tabla vemos que hay algunas conjeturas, o posibles valores de \\(\\theta\\), que tienen probabilidad considerablemente más alta que otra. La notación\n\\[p(D|\\theta)\\] significa: la probabilidad de los datos \\(D\\) dado el valor de \\(\\theta\\). Nótese que esta distribución no es la posterior que describimos arriba, y no es una distribución de probabilidad sobre \\(\\theta\\) (las probabilidades no suman uno). Esta función se llama usualmente verosimilitud de los datos, e incorpora supuestos concretos del proceso generador de los datos.\nUsando reglas de probabilidad (en particular la regla de Bayes), observamos que\n\\[p(\\theta | D) = \\frac{p(D|\\theta)p(\\theta)} { p(D)}.\\] Como \\(p(\\theta|D)\\) debe dar una distribución de probabilidad (suma o integra a 1), entonces \\(p(D)\\) debe ser una constante de normalización para el numerador de la derecha, es decir, basta escribir\n\\[p(\\theta | D) \\propto p(D|\\theta)p(\\theta) \\] Ahora es donde encontramos que tenemos que tener \\(p(\\theta)\\) para poder calcular la cantidad que nos interesa, que es la distribución posterior \\(p(\\theta|D)\\). \\(p(\\theta)\\), la distribución a priori o distribución inicial es simplemente una afirmación de dónde puede estar \\(\\theta\\), antes de observar ningún dato.\nPor el momento, podríamos poner \\(p(\\theta)\\) constante, de manera que es parte de la constante de normalización, y sólo tendríamos que normalizar como sigue:\n\ntheta <- seq(0, 1, length.out = 11)\nprob_post <- tibble(conjetura = theta, probablidad = theta^2 * (1 - theta)^3) |> \n mutate(prob_posterior = probablidad / sum(probablidad)) \nprob_post |> \n kable(col.names = c(\"Conjetura θ\", \"p(D|θ)\",\"p(θ|D)\")) |>\n kable_paper()\n\n\n\n\nConjetura θ\np(D|θ)\np(θ|D)\n\n\n\n\n0.0\n0.00000\n0.0000000\n\n\n0.1\n0.00729\n0.0437444\n\n\n0.2\n0.02048\n0.1228923\n\n\n0.3\n0.03087\n0.1852385\n\n\n0.4\n0.03456\n0.2073807\n\n\n0.5\n0.03125\n0.1875188\n\n\n0.6\n0.02304\n0.1382538\n\n\n0.7\n0.01323\n0.0793879\n\n\n0.8\n0.00512\n0.0307231\n\n\n0.9\n0.00081\n0.0048605\n\n\n1.0\n0.00000\n0.0000000\n\n\n\n\n\n\n\nCon esto, expresamos nuestro conocimiento acerca de \\(\\theta\\), después de observar los datos, con una distribución posterior de probabilidad sobre las posibles conjecturas. Este es el resultado principal de inferencia bayesiana, y es la base para tomar decisiones relativas a \\(\\theta\\).\n\nUsando información adicional\nSupongamos que tenemos información adicional acerca de \\(\\theta\\), por ejemplo, que en un experimento similar anterior alguien tomó una muestra de dos personas, y encontraron dos negativos. Tenemos entonces como creencias inciales:\n\ntheta <- seq(0, 1, length.out = 11)\nprob_priori <- tibble(conjetura = theta) |> \n mutate(prob_priori = (1 - theta) * (1 - theta)) |> \n mutate(prob_priori = prob_priori / sum(prob_priori)) \nprob_priori |>\n kable(col.names = c(\"Conjetura θ\", \"p(θ)\")) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\n\n\n\n\n0.0\n0.2597403\n\n\n0.1\n0.2103896\n\n\n0.2\n0.1662338\n\n\n0.3\n0.1272727\n\n\n0.4\n0.0935065\n\n\n0.5\n0.0649351\n\n\n0.6\n0.0415584\n\n\n0.7\n0.0233766\n\n\n0.8\n0.0103896\n\n\n0.9\n0.0025974\n\n\n1.0\n0.0000000\n\n\n\n\n\n\n\nPor ejemplo, al probabilidad inicial de que \\(\\theta\\) sea muy grande es cercana a cero, pues observamos dos negativos y ningún positivo. Ahora regresamos a considerar nuestra fórmula\n\\[p(\\theta | D) \\propto p(D|\\theta)p(\\theta), \\]\nEn este caso, la apriori o inicial tiene un efecto sobre la posterior. Reconsideramos entonces la posterior de nuestra muestra de 5 personas, y calculamos el producto de \\(P(D|\\theta)\\) por \\(p(\\theta)\\):\n\nprob_post <- prob_priori |> \n mutate(verosimilitud = theta^2 * (1 - theta)^3) |> \n mutate(prod = verosimilitud * prob_priori)\n\nprob_post|> \n kable(col.names = c(\"Conjetura θ\", \"p(θ)\", \"p(D|θ)\",\n \"p(D|θ)p(θ)\")) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\np(D|θ)\np(D|θ)p(θ)\n\n\n\n\n0.0\n0.2597403\n0.00000\n0.0000000\n\n\n0.1\n0.2103896\n0.00729\n0.0015337\n\n\n0.2\n0.1662338\n0.02048\n0.0034045\n\n\n0.3\n0.1272727\n0.03087\n0.0039289\n\n\n0.4\n0.0935065\n0.03456\n0.0032316\n\n\n0.5\n0.0649351\n0.03125\n0.0020292\n\n\n0.6\n0.0415584\n0.02304\n0.0009575\n\n\n0.7\n0.0233766\n0.01323\n0.0003093\n\n\n0.8\n0.0103896\n0.00512\n0.0000532\n\n\n0.9\n0.0025974\n0.00081\n0.0000021\n\n\n1.0\n0.0000000\n0.00000\n0.0000000\n\n\n\n\n\n\n\nY finalmente, normalizamos para encontrar la probabilidad posterior:\n\nprob_post <- prob_post |> \n mutate(prob_posterior = prod / sum(prod))\n\nprob_post|> \n kable(col.names = c(\"Conjetura θ\", \"p(θ)\", \"p(D|θ)\",\n \"p(D|θ)p(θ)\", \"p(θ|D)\"), escape = FALSE) |> kable_paper()\n\n\n\n\nConjetura θ\np(θ)\np(D|θ)\np(D|θ)p(θ)\np(θ|D)\n\n\n\n\n0.0\n0.2597403\n0.00000\n0.0000000\n0.0000000\n\n\n0.1\n0.2103896\n0.00729\n0.0015337\n0.0992712\n\n\n0.2\n0.1662338\n0.02048\n0.0034045\n0.2203539\n\n\n0.3\n0.1272727\n0.03087\n0.0039289\n0.2542983\n\n\n0.4\n0.0935065\n0.03456\n0.0032316\n0.2091640\n\n\n0.5\n0.0649351\n0.03125\n0.0020292\n0.1313412\n\n\n0.6\n0.0415584\n0.02304\n0.0009575\n0.0619745\n\n\n0.7\n0.0233766\n0.01323\n0.0003093\n0.0200177\n\n\n0.8\n0.0103896\n0.00512\n0.0000532\n0.0034430\n\n\n0.9\n0.0025974\n0.00081\n0.0000021\n0.0001362\n\n\n1.0\n0.0000000\n0.00000\n0.0000000\n0.0000000\n\n\n\n\n\n\n\nLa última columna nos da el resultado final de la inferencia bayesiana. Podemos resumir algunas de sus características, por ejemplo:\n\nEs muy poco probable que la seropositividad sea mayor o igual a 0.7\nUn intervalo de 90% de probabilidad para la seropositividad es \\([0.1, 0.5]\\)\n\nLa gráfica de la posterior es:\n\nprob_post |>\n ggplot(aes(x = conjetura, y = prob_posterior)) +\n geom_col() +\n labs(x = \"theta\", y = \"Prob posterior\") \n\n\n\n\nAhora podemos definir, para nuestro ejemplo discretizado, la función que calcula la posterior dados los pasos 1 y 2:\n\ncalcular_posterior <- 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 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\n\nmuestra <- c(1,0,0,1,0)\n\n\ncalcular_posterior(muestra, prob_priori) \n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0.0993 \n 3 0.2 0.220 \n 4 0.3 0.254 \n 5 0.4 0.209 \n 6 0.5 0.131 \n 7 0.6 0.0620 \n 8 0.7 0.0200 \n 9 0.8 0.00344 \n10 0.9 0.000136\n11 1 0 \n\n\nProcedemos ahora a hacer algunas pruebas simples de nuestra función:\n\ncalcular_posterior(rep(0, 50)) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0.996\n 2 0.1 0.004\n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0 \n 6 0.5 0 \n 7 0.6 0 \n 8 0.7 0 \n 9 0.8 0 \n10 0.9 0 \n11 1 0 \n\ncalcular_posterior(rep(1, 50)) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0 \n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0 \n 6 0.5 0 \n 7 0.6 0 \n 8 0.7 0 \n 9 0.8 0.011\n10 0.9 0.989\n11 1 0 \n\ncalcular_posterior(c(rep(0, 100), rep(1, 100))) |> round(3)\n\n# A tibble: 11 × 2\n theta prob_posterior\n <dbl> <dbl>\n 1 0 0 \n 2 0.1 0 \n 3 0.2 0 \n 4 0.3 0 \n 5 0.4 0.023\n 6 0.5 0.966\n 7 0.6 0.01 \n 8 0.7 0 \n 9 0.8 0 \n10 0.9 0 \n11 1 0 \n\n\n\n\nMás verificaciones a priori\nOtra verificación útil que podemos hacer es, una vez que hemos definido nuestro modelo generativo y un modelos estadístico asociado, generar bajo simulación datos que podríamos observar. Esto tiene como fin verificar que nuestro modelo generativo y nuestro modelo estadístico producen datos que están de acuerdo con el conocimiento experto (teoría científica o conocimiento de negocio).\nAsí que simulamos datos del modelo:\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \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 tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\n\nPodemos ver por ejemplo dónde esperamos ver el número de positivos a lo largo de distintas muestras, cuando \\(N=30\\):\n\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos), .groups = \"drop\") |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n\n\n\nObservamos que con nuestros supuestos, hay una probabilidad alta de observar 0 positivos (alrededor de 0.30). Esto se debe en parte a la discretización que hicimos, y que nuestra apriori pone peso considerable en prevalencia igual a cero, lo que quizá no es muy realista, y probablemente deberíamos escoger al menos una discretización más fina.\nTambién, si consideramos los supuestos como correctos, esto puede indicar el riesgo de usar una muestra chica para estimar prevalencia si esta es muy baja: es probable que obtengamos 0 observaciones positivas.\n\n\n\n\n\n\nVerificación predictiva a priori\n\n\n\nCon este tipo de verificaciones podemos detectar las consecuencias de nuestros supuestos (incluyendo la elección de distribuciones a priori), así como otras decisiones de modelado (como la discretización).\nConflictos con el conocimiento del área deben ser explorados para entenderlos y si es necesario corregir nuestros supuestos.\n\n\nEste tipo de verificaciones es muy flexible, y debe adaptarse a los aspectos del conocimiento del área que son importantes para los expertos. Podemos usar todos nuestros recursos analíticos (tablas, resúmenes, gráficas) para producir estos chequeos." }, { "objectID": "02-flujo-basico.html#paso-4-probar-el-proceso-de-estimación", "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 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." + "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 como coinciden o no datos generados con la 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." }, { "objectID": "02-flujo-basico.html#paso-5-analizar-los-datos-y-resumir-resultados.", @@ -88,21 +88,14 @@ "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 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\nMás de verificaciones apriori\nAntes de continuar, sin embargo, veremos cómo se veo el chequeo predictivo a priori que consideramos en la sección de arriba.\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \n function(rep){\n # apriori seleccionada\n theta_sim <- rbeta(1, 1, 3)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\n\nPodemos ver por ejemplo dónde esperamos ver el número de positivos a lo largo de distintas muestras, cuando \\(N=30\\):\n\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos)) |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n`summarise()` has grouped output by 'rep'. You can override using the `.groups`\nargument.\n\n\n\n\n\nEste resultado es consecuencia de nuestros supuestos, antes de ver los datos, y resume que esperamos con mayor probabilidad un número bajo de positivos (en una muestra de N=30), y que es muy poco probable que observemos prevalencias muy altas. Dependiendo de la situación, este puede ser un resultado aceptable.\nUn resultado no aceptable para una enfermedad que sabemos que es relativamente rara (aunque tenemos incertidumbre), por ejemplo, podría ser el siguiente:\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \n function(rep){\n # apriori seleccionada\n theta_sim <- rbeta(1, 30, 3)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos)) |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n`summarise()` has grouped output by 'rep'. You can override using the `.groups`\nargument.\n\n\n\n\n\nEste resultado no es aceptable cuando sabemos que es prácticamente imposible que la mayoría de la población está infectada. Debemos entonces regresar y ajustar nuestros supuestos: el problema en este caso es la elección de la distribución a priori para \\(\\theta\\).\nObservación: la crítica es sobre el conjunto completo de supuestos iniciales que hacemos acerca del problema. Cuando los diagnósticos no son aceptables desde el punto de vista teórico es necesario investigar dónde está el problema. Las distribuciones apriori que usamos, igual que cualquier supuesto, están sujetas a esta crítica. Nótese que esta crítica la estamos haciendo sin ver los datos que esperamos observar: es una crítica de supuestos.\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.4280916\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.0961\n\n\nMuchas veces se presentan intervalos de probabilidad posterior, por ejemplo, podríamos reportar que con 90% de probabilidad la prevalencia está en el siguiente intervalo:\n\nsims_post |> \n summarise(inf = quantile(theta, 0.05),\n sup = quantile(theta, 0.95)) |> \n mutate(inf = round(inf, 2),\n sup = round(sup, 2))\n\n# A tibble: 1 × 2\n inf sup\n <dbl> <dbl>\n1 0.16 0.73\n\n\nObservación: No hay un intervalo mágico que debe reportarse (por ejemplo 95% de probabilida es una costumbre o superstición). Hay varias maneras de construir intervalos de probabilidad. Dejaremos esta discusión para más adelante.\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)" + "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 ahora. 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}(1-\\theta)^{b-1}\\]\nentonces la posterior, por la fórmula de Bayes, es:\n\\[p(\\theta|D) \\propto \\theta^{N_+ +a -1 }(1-\\theta)^{N_{-}+b-1}\\] que también es de la familia beta, pero con parámetros \\((N_{+} +a, N_{-} +b)\\).\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, \n densidad = dbeta(theta_seq, n_pos + 1, n_neg + 1)))) |> \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, n_pos + 1, n_neg + 3)))) |> \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\nMás de verificaciones apriori\nAntes de continuar, sin embargo, veremos cómo se veo el chequeo predictivo a priori que consideramos en la sección de arriba.\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \n function(rep){\n # apriori seleccionada\n theta_sim <- rbeta(1, 1, 3)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\n\nPodemos ver por ejemplo dónde esperamos ver el número de positivos a lo largo de distintas muestras, cuando \\(N=30\\):\n\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos)) |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n`summarise()` has grouped output by 'rep'. You can override using the `.groups`\nargument.\n\n\n\n\n\nEste resultado es consecuencia de nuestros supuestos, antes de ver los datos, y resume que esperamos con mayor probabilidad un número bajo de positivos (en una muestra de N=30), y que es muy poco probable que observemos prevalencias muy altas. Dependiendo de la situación, este puede ser un resultado aceptable.\nUn resultado no aceptable para una enfermedad que sabemos que es relativamente rara (aunque tenemos incertidumbre), por ejemplo, podría ser el siguiente:\n\nset.seed(231)\nsimulacion_datos_tbl <- map_df(1:500, \n function(rep){\n # apriori seleccionada\n theta_sim <- rbeta(1, 30, 3)\n datos_sim <- sim_pos_neg(theta = theta_sim, N = 30)\n tibble(rep = rep, theta_sim = theta_sim, datos_sim)\n })\nsimulacion_datos_tbl |> \n group_by(rep, theta_sim) |> \n summarise(Npos = sum(Pos)) |> \n ggplot(aes(x = Npos)) +\n geom_bar() +\n labs(x = \"Número de positivos\", y = \"Frecuencia (muestras)\") \n\n`summarise()` has grouped output by 'rep'. You can override using the `.groups`\nargument.\n\n\n\n\n\nEste resultado no es aceptable cuando sabemos que es prácticamente imposible que la mayoría de la población está infectada. Debemos entonces regresar y ajustar nuestros supuestos: el problema en este caso es la elección de la distribución a priori para \\(\\theta\\).\nObservación: la crítica es sobre el conjunto completo de supuestos iniciales que hacemos acerca del problema. Cuando los diagnósticos no son aceptables desde el punto de vista teórico es necesario investigar dónde está el problema. Las distribuciones apriori que usamos, igual que cualquier supuesto, están sujetas a esta crítica. Nótese que esta crítica la estamos haciendo sin ver los datos que esperamos observar: es una crítica de supuestos.\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.4280916\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.0961\n\n\nMuchas veces se presentan intervalos de probabilidad posterior, por ejemplo, podríamos reportar que con 90% de probabilidad la prevalencia está en el siguiente intervalo:\n\nsims_post |> \n summarise(inf = quantile(theta, 0.05),\n sup = quantile(theta, 0.95)) |> \n mutate(inf = round(inf, 2),\n sup = round(sup, 2))\n\n# A tibble: 1 × 2\n inf sup\n <dbl> <dbl>\n1 0.16 0.73\n\n\nObservación: No hay un intervalo mágico que debe reportarse (por ejemplo 95% de probabilidad es una costumbre o superstición). Hay varias maneras de construir intervalos de probabilidad. Dejaremos esta discusión para más adelante.\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\\) en general 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\\), usaríamos \\(f(\\theta) = \\theta\\). 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 un número grande \\(M\\) de valores de la posterior (aquí usaremos \\(M=10000\\))\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 })\nsimular_posterior <- function(muestra, n){\n tibble(theta = \n rbeta(n, sum(muestra) + 1, \n length(muestra) - sum(muestra) + 3))\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", "href": "02-flujo-basico.html#observaciones-1", "title": "2  Flujo de trabajo básico: motivación", "section": "2.8 Observaciones", - "text": "2.8 Observaciones\nEl proceso de arriba lo refinaremos considerablemente en el resto del curso.\n\nEn primer lugar, los modelos generativos serán más complicados, y estarán basados en teoría más compleja (que expresamos con diagramas causales)\nUsaremos más herramientas y componentes para construir modelos estadísticos apropiados, ya sea que construyamos un modelo completo para todo el proceso de generación de datos, o que usemos modelos estándar como regresión para aproximar respuestas, cuando es apropiado\nRefinaremos el proceso de checar que el cómputo (checar Monte Carlo) y la inferencia (verificación apriori) es correcta bajo nuestros supuestos.\nFinalmente, veremos qué hacer después de hacer la estimación y que los puntos de arriba están resueltos, para tener confianza en nuestras conclusiones." - }, - { - "objectID": "02-flujo-basico.html#ventajas-y-desventajas-de-métodos-bayesianos", - "href": "02-flujo-basico.html#ventajas-y-desventajas-de-métodos-bayesianos", - "title": "2  Flujo de trabajo básico: motivación", - "section": "2.9 Ventajas y desventajas de métodos bayesianos", - "text": "2.9 Ventajas y desventajas de métodos bayesianos" + "text": "2.8 Observaciones\nEl proceso de arriba lo refinaremos considerablemente en el resto del curso.\n\nEn primer lugar, los modelos generativos serán más complicados, y estarán basados en teoría más compleja (que expresamos con diagramas causales)\nUsaremos más herramientas y componentes para construir modelos estadísticos apropiados, ya sea que construyamos un modelo completo para todo el proceso de generación de datos, o que usemos modelos estándar como regresión para aproximar respuestas, cuando es apropiado\nRefinaremos el proceso de checar que el cómputo (checar Monte Carlo) y la inferencia (verificación apriori) es correcta bajo nuestros supuestos.\nFinalmente, veremos qué hacer después de hacer la estimación y que los puntos de arriba están resueltos, para tener confianza en nuestras conclusiones.\n\n\n2.8.1 Resumen\nAquí juntamos algunas observaciones que se derivan de lo que hemos visto (flujo de trabajo y estimación bayesiana):\n\nTodo nuestro trabajo está fundamentado en entender qué es lo que queremos estimar dentro de un modelo generativo. Los diagramas causales nos ayudan a conectar el problema de interés con nuestros modelos, a construir modelos generativos y hacer explícitos nuestros supuestos.\nEl proceso de estimación siempre es el mismo: nuestro estimador es la distribución posterior, que se construye a partir de la verosimilitud y la apriori (modelo generativo). Nuestro estimador es la posterior de las cantidades de interés, que pueden resumirse de distintas maneras. Cualquier cálculo derivado de otras cantidades de interés debe considerar toda la posterior (no solo la media o la moda, etc. posterior).\nNuestro proceso incluye los chequeos predictivos a priori (basados en simulación de datos). Esto son cruciales para detectar problemas en nuestros supuestos (vs teoría) y que nuestro proceso sea internamente consistente. Esto también es una verificación de la información a priori.\nGeneralmente es más conveniente y práctico hacer simulaciones que calcular analíticamente la posterior o sus integrales." }, { "objectID": "02-flujo-basico-2.html#prevalencia-con-error-conocido",