Este es un curso de modelación bayesiana aplicada, que se concentra en plantear y resolver problemas aplicados usando estadística. Para hacer esto necesitamos entender tres componentes:
El contexto del problema, o conocimiento del área: esto es conocimiento científico o de negocios y cómo es que observamos o medimos el fenómeno que nos interesa.
En vista del contexto del problema (1), ¿cómo construir modelos estadísticos para responder preguntas de interés?
Práctica de construcción y desarrollo de modelos, junto con procedimientos para criticarlos y probarlos.
Este proceso aplica tanto a estadística bayesiana como frecuentista, aunque en este curso, por su flexibilidad y unidad, consideraremos el enfoque bayesiano. Para la parte 1 y 2, consideraremos modelos gráficos causales, que expresan nuestro conocimiento acerca del problema de interés. Para la parte 3 propondremos un flujo de trabajo bayesiano que nos permita probar y entender el funcionamiento de nuestros modelos. Finalmente, para 2 y 3 propondremos distintas estrategias de modelación, como son modelos jerárquicos, regresión, como trabajar con datos faltantes entre otros.
En primer lugar, observamos (McElreath (2020)):
Las 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.
Las 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.).
Consideremos un ejemplo simple para ilustrar este primer principio:
Este 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.
La tabla original tiene 700 renglones (cada renglón es un paciente)
calculos <- read_csv("../datos/kidney_stone_data.csv")
names(calculos) <- c("tratamiento", "tamaño", "éxito")
calculos <- calculos |>
mutate(tamaño = ifelse(tamaño == "large", "grandes", "chicos")) |>
mutate(resultado = ifelse(éxito == 1, "mejora", "sin_mejora")) |>
select(tratamiento, tamaño, resultado)
nrow(calculos)[1] 700y se ve como sigue (muestreamos algunos renglones):
calculos |>
sample_n(10) |> kable() |>
kable_paper(full_width = FALSE)| tratamiento | tamaño | resultado |
|---|---|---|
| B | chicos | sin_mejora |
| B | chicos | mejora |
| B | chicos | mejora |
| B | chicos | mejora |
| A | grandes | mejora |
| A | grandes | mejora |
| A | grandes | mejora |
| A | grandes | mejora |
| A | grandes | mejora |
| A | chicos | mejora |
Aunque estos datos contienen información de 700 pacientes, los datos pueden resumirse sin pérdida de información contando como sigue:
calculos_agregada <- calculos |>
group_by(tratamiento, tamaño, resultado) |>
count()
calculos_agregada |> kable() |>
kable_paper(full_width = FALSE)| tratamiento | tamaño | resultado | n |
|---|---|---|---|
| A | chicos | mejora | 81 |
| A | chicos | sin_mejora | 6 |
| A | grandes | mejora | 192 |
| A | grandes | sin_mejora | 71 |
| B | chicos | mejora | 234 |
| B | chicos | sin_mejora | 36 |
| B | grandes | mejora | 55 |
| B | grandes | sin_mejora | 25 |
Como en este caso nos interesa principalmente la tasa de éxito de cada tratamiento, podemos mejorar mostrando como sigue:
calculos_agregada |> pivot_wider(names_from = resultado, values_from = n) |>
mutate(total = mejora + sin_mejora) |>
mutate(prop_mejora = round(mejora / total, 2)) |>
select(tratamiento, tamaño, total, prop_mejora) |>
arrange(tamaño) |>
kable() |>
kable_paper(full_width = FALSE)| tratamiento | tamaño | total | prop_mejora |
|---|---|---|---|
| A | chicos | 87 | 0.93 |
| B | chicos | 270 | 0.87 |
| A | grandes | 263 | 0.73 |
| B | grandes | 80 | 0.69 |
Esta tabla descriptiva es una reescritura de los datos, y no hemos resumido nada todavía. Pero es apropiada para empezar a contestar la pregunta:
Supongamos que otro analista decide comparar los pacientes que recibieron cada tratamiento, ignorando la variable de tamaño:
calculos |> group_by(tratamiento) |>
summarise(prop_mejora = mean(resultado == "mejora") |> round(2)) |>
kable() |>
kable_paper(full_width = FALSE)| tratamiento | prop_mejora |
|---|---|
| A | 0.78 |
| B | 0.83 |
y 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.
Podemos 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:
calculos |> group_by(tratamiento, tamaño) |> count() |>
kable() |>
kable_paper(full_width = FALSE)| tratamiento | tamaño | n |
|---|---|---|
| A | chicos | 87 |
| A | grandes | 263 |
| B | chicos | 270 |
| B | grandes | 80 |
Nuestra 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.
Los resúmenes descriptivos acompañados de hipótesis causales acerca del proceso generador de datos, nos guía hacia descripciones interpretables de los datos.
Las explicaciones no son tan simples y, otra vez, interviene el comportamiento de doctores, tratamientos, y distintos tipos de padecimientos.
Podemos codificar la información causal con un diagrama:
grViz("
digraph {
graph [ranksep = 0.2]
node [shape=plaintext]
T
M
C
edge [minlen = 3]
T -> M
C -> T
C -> M
{ rank = same; M; T }
}
", 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).
Contrastemos 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:
corazon <- calculos |>
select(tratamiento, presión = tamaño, resultado) |>
mutate(presión = ifelse(presión == "grandes", "alta", "baja"))
corazon_agregada <- corazon |>
group_by(tratamiento, presión, resultado) |>
count()
corazon_agregada |> pivot_wider(names_from = resultado, values_from = n) |>
mutate(total = mejora + sin_mejora) |>
mutate(prop_mejora = round(mejora / total, 2)) |>
select(tratamiento, presión, total, prop_mejora) |>
arrange(presión) |>
kable() |>
kable_paper(full_width = FALSE)| tratamiento | presión | total | prop_mejora |
|---|---|---|---|
| A | alta | 263 | 0.73 |
| B | alta | 80 | 0.69 |
| A | baja | 87 | 0.93 |
| B | baja | 270 | 0.87 |
corazon |> group_by(tratamiento) |>
summarise(prop_mejora = mean(resultado == "mejora") |> round(2)) |>
kable() |>
kable_paper(full_width = FALSE)| tratamiento | prop_mejora |
|---|---|
| A | 0.78 |
| B | 0.83 |
¿Cuál creemos que es el mejor tratamiento en este caso? ¿Deberíamos usar la tabla agregada o la desagregada por presión?
Nuestros supuestos causales podemos mostrarlos con el siguiente diagrama:
grViz("
digraph {
graph [ranksep = 0.2]
node [shape=plaintext]
P
T
M
edge [minlen = 3]
T -> P
P -> M
T -> M
{ rank = same; M; T}
}
", 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.
Los diagramas de arriba se llaman DAGs (Gráficas dirigidas acíclicas), y no son generadas por datos observados, sino que codifican conocimiento acerca del fenómenos y los datos observados. Nos ayudan a (McElreath (2020)):
Los DAGs se construyen con causas, e implican asociaciones observables, pero no se construyen con asociaciones simplemente. El pensamiento causal es útil siempre que queremos responder preguntas acerca de un fenómeno de interés. En particular nos asisten en :
Una ejemplo al que regresaremos más adelante es el siguiente: en 1973 se recolectaron datos agregados de solicitantes para estudiar en Berkeley para los 6 departamentos más grandes, clasificados por sexo del solicitante y si fue admitido o no. Los resultados se muestran a continuación:
data("UCBAdmissions")
adm_original <- UCBAdmissions |> as_tibble() |>
pivot_wider(names_from = Admit, values_from = n)
adm_original |> knitr::kable() |>
kable_paper(full_width = FALSE)| Gender | Dept | Admitted | Rejected |
|---|---|---|---|
| Male | A | 512 | 313 |
| Female | A | 89 | 19 |
| Male | B | 353 | 207 |
| Female | B | 17 | 8 |
| Male | C | 120 | 205 |
| Female | C | 202 | 391 |
| Male | D | 138 | 279 |
| Female | D | 131 | 244 |
| Male | E | 53 | 138 |
| Female | E | 94 | 299 |
| Male | F | 22 | 351 |
| Female | F | 24 | 317 |
y las proporciones de admisión por sexo y departamente son las siguientes:
adm_tbl <- adm_original |>
mutate(prop_adm = round(Admitted / (Admitted + Rejected), 2), total = Admitted + Rejected) |>
select(Gender, Dept, prop_adm, total) |>
pivot_wider(names_from = Gender, values_from = prop_adm:total)
adm_tbl |> knitr::kable() |>
kable_paper(full_width = FALSE)| Dept | prop_adm_Male | prop_adm_Female | total_Male | total_Female |
|---|---|---|---|---|
| A | 0.62 | 0.82 | 825 | 108 |
| B | 0.63 | 0.68 | 560 | 25 |
| C | 0.37 | 0.34 | 325 | 593 |
| D | 0.33 | 0.35 | 417 | 375 |
| E | 0.28 | 0.24 | 191 | 393 |
| F | 0.06 | 0.07 | 373 | 341 |
Complementamos con las tasas de aceptación a total por género, y tasas de aceptación por departamento:
adm_original |> group_by(Gender) |>
summarise(Admitted = sum(Admitted), Rejected = sum(Rejected)) |>
mutate(prop_adm = round(Admitted / (Admitted + Rejected),2)) |>
kable() |>
kable_paper(full_width = FALSE)| Gender | Admitted | Rejected | prop_adm |
|---|---|---|---|
| Female | 557 | 1278 | 0.30 |
| Male | 1198 | 1493 | 0.45 |
La pregunta que queremos hacer es: ¿existe discriminación por sexo en la selección de candidatos? Examinando las tablas no está clara cuál es la respuesta.
adm_original |> group_by(Dept) |>
summarise(Admitted = sum(Admitted), Rejected = sum(Rejected)) |>
mutate(prop_adm = round(Admitted / (Admitted + Rejected),2)) |>
kable() |>
kable_paper(full_width = FALSE)| Dept | Admitted | Rejected | prop_adm |
|---|---|---|---|
| A | 601 | 332 | 0.64 |
| B | 370 | 215 | 0.63 |
| C | 322 | 596 | 0.35 |
| D | 269 | 523 | 0.34 |
| E | 147 | 437 | 0.25 |
| F | 46 | 668 | 0.06 |
Discutiremos este ejemplo con más detalle más adelante. La interpretación debe ser hecha con cuidado, y debemos establecer claramente los supuestos que fundamentan nuestra decisión de mostrar cada tabla y de qué forma mostrarlas.
En muchos cursos introductorios de estadística se muestran distintos tipos de procedimientos, que aplican según el tipo de datos (por ejemplo, categóricos o numéricos, pareados, no pareados, etc), generalmente con el propósito de evaluar evidencia en contra de una hipótesis nula. Por ejemplo, de McElreath (2020):
Este enfoque puede ser confuso en un principio (¿cómo se relacionan todos estos procedimientos?), y también restringir nuestra capacidad para analizar datos: ¿qué hacemos cuando no se cumplen los supuestos de un procedimiento? Adicionalmente si no tenemos mucha experiencia, la manera en que fallan estas herramientas puede ser poco intuitiva y difícil de descubrir.
Y aunque son herramientas poderosas, no sustituyen el pensamiento científico o de proceso de negocios. Estas herramientas no generan hallazgos si no están acompañados de pensamiento causal.
Buscamos entonces:
Iremos refinando nuestro poco a poco, conforme veamos distintas herramientas y problemas. El más básico es el siguiente (McElreath (2020)):
Este proceso no es exclusivo de los modelos bayesianos, pero quizá es más natural, como veremos, cuando adoptamos el punto de vista bayesiano. Su propósito es múltiple: verificar que nuestros modelos están estimando las cantidades que realmente nos interesan, según nuestros supuestos, verificar los programas y cómputos con los que se obtienen resultados, y checar la adecuación del modelo a datos reales, cuestionando supuestos teóricos y supuestos de modelación.
Finalmente, quisiéramos llegar a un proceso como el que se describe en Towards a Principled Bayesian Workflow, e incorporar el que se detalla en Gelman et al. (2020):
Cualquier proceso de análisis de datos se beneficia de muchos aspectos de ingenería de software. Parte de la profesionalización del análisis de datos que observamos en ciencia de datos es utilizar las herramientas reconocidas para resolver problemas de desarrollo y calidad de código, así como su documentación.
Análisis como software: Una parte de este proceso está relacionado con la reproducibilidad y documentación del trabajo, y su objetivo es evitar errores de programación y de organización (esta parte hablaremos menos: es necesario seguir los estándares de la industria para obtener resultados más confiables).
Otra parte es el proceso con el cual construimos y contrastamos modelos para contestar preguntas, verificamos los modelos y sus respuestas y checamos resultados de cómputos.