Código
library(tidyverse)
library(knitr)
library(kableExtra)
library(ggplot2)Analítica Avanzada de Datos
Challenge #1
Introducción
En esta actividad se analizarán datos provenientes de una tarea repetitiva en niños y adultos diagnosticados y no diagnosticados en Trastorno de Déficit de Atención con Hiperactividad (ADHD en inglés), residentes en Barranquilla.
Los datos se encuentran aquí.
Las variables de interés son sexo, edad, adhdstatus y las variables del seguimiento trait_A1, trait_A2, \(...\), trait_A50. Estas últimas representan el estado de finalización de la tarea. La tarea se cuantifica en una escala de 1 a 5. Por simplicidad, valores diferentes a 5 se consideran no exitosos.
Fecha de entrega: Sábado, 25 de Mayo de 2024.
Preguntas
P1
Suponga que la respuesta de un individuo puede resumirse en el vector 01000111111101100100001110100111100110000111001110. Construya una función que permita calcular la entropía de Shannon a partir de un string de 0s y 1s. Calcule la entropía \(H\) para cada individuo. Estudie la distribución de \(H_1, H_2,\ldots, H_n\) por sexo y dignóstico de ADHD. Qué observa?
Respuesta.
Cargamos las librerias necesarias
Inicialmente cargamos los datos
Código
datos <- read.table("/Users/josejavier-yenifer/Documents/Analitica avanzada de Datos/taskA.txt", header = TRUE, sep = "\t")
knitr::kable(head(datos))| uid | sexo | edad | adhdstatus | ina | hyp | combined | trait_A1 | trait_A2 | trait_A3 | trait_A4 | trait_A5 | trait_A6 | trait_A7 | trait_A8 | trait_A9 | trait_A10 | trait_A11 | trait_A12 | trait_A13 | trait_A14 | trait_A15 | trait_A16 | trait_A17 | trait_A18 | trait_A19 | trait_A20 | trait_A21 | trait_A22 | trait_A23 | trait_A24 | trait_A25 | trait_A26 | trait_A27 | trait_A28 | trait_A29 | trait_A30 | trait_A31 | trait_A32 | trait_A33 | trait_A34 | trait_A35 | trait_A36 | trait_A37 | trait_A38 | trait_A39 | trait_A40 | trait_A41 | trait_A42 | trait_A43 | trait_A44 | trait_A45 | trait_A46 | trait_A47 | trait_A48 | trait_A49 | trait_A50 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | M | 35 | 1 | 0 | 0 | 1 | 1 | 5 | 4 | 4 | 4 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 4 | 5 | 5 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 5 | 5 | 5 | 4 | 5 | 4 | 4 | 5 | 5 | 5 | 5 | 4 | 4 | 5 | 5 | 4 | 4 | 4 | 4 | 5 | 5 | 5 | 4 | 4 | 5 | 5 | 5 | 4 |
| 2 | F | 39 | 0 | 0 | 0 | 0 | 5 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 4 | 4 | 4 | 4 | 4 | 4 | 5 | 4 | 4 | 5 | 4 | 5 | 4 | 5 | 5 | 4 | 5 | 5 | 4 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 4 | 4 | 4 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 5 | 5 | 4 | 4 | 4 | 5 |
| 3 | F | 7 | 0 | 0 | 0 | 0 | 1 | 5 | 4 | 4 | 5 | 5 | 5 | 5 | 4 | 3 | 3 | 5 | 5 | 5 | 5 | 5 | 4 | 5 | 5 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 4 | 5 | 3 | 3 | 5 | 4 | 1 | 3 | 5 | 4 | 5 | 1 | 5 | 4 | 4 | 4 | 4 | 4 | 5 | 1 | 5 | 4 | 4 | 5 |
| 8 | M | 38 | 1 | 1 | 0 | 0 | 5 | 5 | 5 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 5 | 4 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 5 | 4 | 5 | 1 | 5 | 5 | 5 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 4 | 4 | 4 |
| 9 | F | 42 | 1 | 1 | 0 | 0 | 4 | 5 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 4 | 5 | 4 | 5 | 5 | 4 | 4 | 4 | 4 | 5 | 4 | 4 | 5 | 4 | 5 | 4 | 5 | 4 | 4 | 4 | 5 | 5 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 5 | 4 | 4 | 4 | 5 | 5 | 4 | 5 | 4 | 4 |
| 10 | M | 12 | 0 | 0 | 0 | 0 | 4 | 3 | 5 | 5 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 5 | 5 | 5 | 4 | 5 | 5 | 5 | 5 | 4 | 5 | 4 | 5 | 3 | 4 | 2 | 4 | 5 | 5 | 5 | 5 | 4 | 4 | 4 | 5 | 4 | 4 | 5 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 5 | 2 | 4 | 5 |
Convertirmos las columnas de seguimiento en un vector binario considerando respuestas diferente a 5 como 0s
Código
datos <- datos %>%
mutate(across(starts_with("trait_A"), ~ as.integer(. == 5)))Código
datos <- datos %>%
rowwise() %>%
mutate(Respuestas = paste(across(starts_with("trait_A")), collapse = "")) %>%
ungroup() Código
knitr::kable(head(datos %>% select(uid, Respuestas)))| uid | Respuestas |
|---|---|
| 1 | 01000111111101100100001110100111100110000111001110 |
| 2 | 11110101000000100101011011011101010000010000110001 |
| 3 | 01001111000111110110010000010010001010100000101001 |
| 8 | 11111101011011111111100100000000010101111110101000 |
| 9 | 01111010101011000010010101000110000000001000110100 |
| 10 | 00110010000111011110101000011110001001000000001001 |
generamos una funcion para calcular la entropía de Shannon \[ \text{H}(X) = -\sum_{x \in \mathcal{X}} p(x) \log_b p(x) \]
Código
calcular_entropia_shannon <- function(string_binario) {
vector_binario <- as.numeric(unlist(strsplit(string_binario, "")))
frecuencias <- table(vector_binario)
probabilidades <- frecuencias / sum(frecuencias)
entropia <- -sum(probabilidades * log2(probabilidades))
return(entropia)
}Se calcula la entropía para cada persona
Código
datos <- datos %>%
mutate(entropia = sapply(Respuestas, calcular_entropia_shannon))Código
knitr::kable(head(datos %>% select(uid, Respuestas,entropia)))| uid | Respuestas | entropia |
|---|---|---|
| 1 | 01000111111101100100001110100111100110000111001110 | 0.9953784 |
| 2 | 11110101000000100101011011011101010000010000110001 | 0.9895875 |
| 3 | 01001111000111110110010000010010001010100000101001 | 0.9814539 |
| 8 | 11111101011011111111100100000000010101111110101000 | 0.9814539 |
| 9 | 01111010101011000010010101000110000000001000110100 | 0.9580420 |
| 10 | 00110010000111011110101000011110001001000000001001 | 0.9709506 |
Generamos un resumen de distribución de la entropía por sexo y diagnóstico de ADHD
Código
datos %>%
group_by(sexo, adhdstatus) %>%
summarise(
media_entropia = mean(entropia),
sd_entropia = sd(entropia),
.groups = 'drop'
) %>%
arrange(desc(media_entropia)) %>%
kable(format = "html", caption = "Distribución de la entropía por sexo y diagnóstico de ADHD") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| sexo | adhdstatus | media_entropia | sd_entropia |
|---|---|---|---|
| F | 1 | 0.9835175 | 0.0204985 |
| F | 0 | 0.9822861 | 0.0215635 |
| M | 1 | 0.9809883 | 0.0249778 |
| M | 0 | 0.9789361 | 0.0356549 |
Código
ggplot(datos, aes(x = interaction(sexo, adhdstatus), y = entropia, fill = sexo)) +
geom_violin(trim = TRUE) +
scale_fill_manual(values = c("pink", "blue")) +
labs(
title = "Distribución de la Entropía por Sexo y Diagnóstico de ADHD",
x = "Grupo",
y = "Entropía",
fill = "Sexo"
) +
theme_minimal() +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5)
)
Los datos revelan diferencias significativas en la entropía de respuesta en tareas repetitivas entre individuos con y sin diagnóstico de ADHD, siendo la entropía media ligeramente superior en aquellos diagnosticados, lo que sugiere una mayor imprevisibilidad en su comportamiento o capacidad de atención.
La variabilidad asociada al ADHD parece afectar más a los hombres que a las mujeres, evidenciado por una mayor diferencia en la entropía media entre estos dos grupos. Además, la desviación estándar de la entropía es más alta en personas sin diagnóstico de ADHD, con una variabilidad especialmente marcada en los hombres, indicando una dispersión más amplia en sus respuestas a las tareas.
Esto podría reflejar una heterogeneidad en la población general respecto a cómo se manejan estas tareas repetitivas.
P2
Conside el vector anterior. Construya una función que calcule el número de cambios de estado. En este caso, una matriz de transición podría ayudar.
Respuesta.
creamos una función en R que tome un vector binario y calcule el número de transiciones de 0 a 1 y de 1 a 0.
Código
calcular_transiciones <- function(string_binario) {
vector_binario <- as.numeric(unlist(strsplit(string_binario, "")))
matriz_transicion <- matrix(0, nrow = 2, ncol = 2, dimnames = list(c("0", "1"), c("0", "1")))
for (i in 2:length(vector_binario)) {
estado_anterior <- vector_binario[i - 1] + 1 # +1 porque los índices en R comienzan en 1
estado_actual <- vector_binario[i] + 1
matriz_transicion[estado_anterior, estado_actual] <- matriz_transicion[estado_anterior, estado_actual] + 1
}
return(matriz_transicion)
}A modo de Ejemplo realizamos la matriz de transición para el primer sujeto (uid=1)
Código
obs_1 <- datos %>% filter(uid == 1) %>% pull(Respuestas)
mat_trans_1 <- calcular_transiciones(obs_1)
knitr::kable(mat_trans_1 )| 0 | 1 | |
|---|---|---|
| 0 | 12 | 10 |
| 1 | 10 | 17 |
P3
Construya la matriz de transición para el individuo 83. Calcule \[P(s+1 = j | s = i)\] donde \(s\) es el intento e $i,j=\{$0,1\(\}\). Defina esta matriz de probabilidades como \(\mathbf{P}^{(1)}\). Tenga en cuenta que \(\mathbf{P}^{(1)}\) es de dimensión \(2\times 2\).
Respuesta.
Definimos:
- \(P^{(1)}_{00}\): Probabilidad de no éxito seguido por otro no éxito.
- \(P^{(1)}_{01}\): Probabilidad de no éxito seguido por éxito.
- \(P^{(1)}_{10}\): Probabilidad de éxito seguido por no éxito.
- \(P^{(1)}_{11}\): Probabilidad de éxito seguido por otro éxito.
Esta matriz será útil para entender la dinámica de comportamiento del individuo 83 en tareas repetitivas y puede ayudar en la planificación de intervenciones dirigidas para mejorar la consistencia en el rendimiento, especialmente en contextos relacionados con ADHD.
Primero extraemo los resultados para la observacion numero 83
Código
obs_83 <- datos %>% filter(uid == 83) %>% pull(Respuestas)
obs_83[1] "01100001011001010111101000101011000100010010010000"ahora aplicamos la funcion para calcular la matriz de transición
Código
mat_trans_83 <- calcular_transiciones(obs_83)
kable(mat_trans_83, caption = "Matriz de transición individuo 83") %>%
kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped", "hover", "condensed")) %>%
row_spec(0, bold = TRUE) %>%
column_spec(1, bold = TRUE)| 0 | 1 | |
|---|---|---|
| 0 | 15 | 14 |
| 1 | 14 | 6 |
Creamos una función para calcular la matriz de probabilidades de transición.
Código
calcular_probabilidades <- function(matriz_transicion) {
sumas_filas <- rowSums(matriz_transicion)
sumas_filas[sumas_filas == 0] <- 1 # Evitar división por cero
matriz_probabilidades <- sweep(matriz_transicion, 1, sumas_filas, FUN="/")
return(matriz_probabilidades)
}La matriz de probabilidades resultante, \(P^{(1)}\), es:
\[ P^{(1)} = \begin{bmatrix} \frac{15}{29} & \frac{14}{29} \\ \frac{14}{20} & \frac{6}{20} \end{bmatrix} \] Calculamos la matriz de probabilidades para el individuo 83
Código
mat_prob_83 <- calcular_probabilidades(mat_trans_83)
kable(mat_prob_83, caption = "Matriz de Probabilidades individuo 83") %>%
kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped", "hover", "condensed")) %>%
row_spec(0, bold = TRUE) %>%
column_spec(1, bold = TRUE)| 0 | 1 | |
|---|---|---|
| 0 | 0.5172414 | 0.4827586 |
| 1 | 0.7000000 | 0.3000000 |
La matriz de transición del individuo 83 muestra una probabilidad de 51.72% de no completar con éxito tareas consecutivas y una casi igual probabilidad de 48.28% de pasar de no completar a completar la tarea, lo que indica una capacidad significativa para mejorar. Una vez que se logra el éxito, hay una alta probabilidad de 70% de mantener ese éxito en el siguiente intento, aunque también existe una probabilidad del 30% de revertir al fracaso.
P4
Construya la matriz de transición en el paso \(k\) utilizando la propiedad de Markov para matrices de transición, esto es
\[ \mathbf{P}^{(k)}_M = \underset{k\text{-veces}}{\underbrace{\mathbf{P}^{(1)}\times \mathbf{P}^{(1)}\times\mathbf{P}^{(1)}\times\cdots\times\mathbf{P}^{(1)}}} \]
En la expresión anterior, \(M\) simboliza el hecho de que está utilizando la propiedad de Markov para matrices de transición. Además, el símbolo \(\times\) se refiere a multiplicación matricial y no a potencias.
Respuesta. Generalizamos lo hecho anteriormente con el fin de que se pueda tomar cualquier sujeto de la misma base de dato o una nueva con la misma estructura.
Código
calcular_P1 <- function(vector_binario) {
matriz_transicion <- matrix(0, nrow = 2, ncol = 2, dimnames = list(c("0", "1"), c("0", "1")))
for (i in 1:(length(vector_binario) - 1)) {
estado_anterior <- vector_binario[i] + 1
estado_siguiente <- vector_binario[i + 1] + 1
matriz_transicion[estado_anterior, estado_siguiente] <- matriz_transicion[estado_anterior, estado_siguiente] + 1
}
sumas_filas <- rowSums(matriz_transicion)
sumas_filas[sumas_filas == 0] <- 1 # Evitar división por cero
matriz_probabilidades <- sweep(matriz_transicion, 1, sumas_filas, FUN = "/")
return(matriz_probabilidades)
}Código
calcular_Pk <- function(uid, k, datos) {
vector_binario <- as.integer(unlist(datos[datos$uid == uid, grep("^trait_A", names(datos))]))
P1 <- calcular_P1(vector_binario)
if (k > 1) {
Pk <- P1
for (j in 2:k) {
Pk <- Pk %*% P1
}
} else {
Pk <- P1
}
return(Pk)
}Probamos con el sujeto número 83 y el paso k=5, de la misma base de datos
Código
P5 <- calcular_Pk(uid = 83, k = 5, datos = datos)
print(P5) 0 1
0 0.5917535 0.4082465
1 0.5919574 0.4080426P5
Calcule \(\mathbf{P}^{(k)}\) como \(P(s+k = j | s = i)\) para \(k=\{1,2,3,4\}\). Use una prueba \(\chi^2\) para tablas de contingencia y compare estos resultados con \(\mathbf{P}^{(k)}_M\). Concluya.
Respuesta. Tomaremos la misma matriz de probabilidad del sujeto N° 83, y a traves de una funcion calculamos las iteraciones para los diferentes K
Código
P1 <- mat_prob_83
# Calcular P^(k) para k = 1, 2, 3, 4
lista_Pk <- lapply(1:4, function(k) {
Pk <- P1
for (j in 2:k) {
Pk <- Pk %*% P1
}
return(Pk)
})luego, vamos a calcular, a traves de una función, $ P^{(k)} $ empíricamente para $ k = 1, 2, 3, 4 $:
Código
calcular_transiciones_empiricas <- function(vector_binario, k) {
n <- length(vector_binario)
matriz <- matrix(0, nrow = 2, ncol = 2, dimnames = list(c("0", "1"), c("0", "1")))
for (i in 1:(n - k)) {
fila <- vector_binario[i] + 1
columna <- vector_binario[i + k] + 1
matriz[fila, columna] <- matriz[fila, columna] + 1
}
suma_filas <- rowSums(matriz)
suma_filas[suma_filas == 0] <- 1
matriz <- sweep(matriz, 1, suma_filas, FUN = "/")
return(matriz)
}Se extrae el vector binario para un individuo (seguimos con el sujeto N° 83)
Código
vector_binario <- as.integer(unlist(datos[datos$uid == 83, grep("^trait_A", names(datos))]))Calculamos y almacenano matrices empíricas para k = 1, 2, 3, 4
Código
lista_Pk_empiricas <- lapply(1:4, function(k) calcular_transiciones_empiricas(vector_binario, k))Para comparar las matrices observacionales con las teóricas $ P^{(k)} $ calculadas anteriormente, usaremos la prueba $ ^2 $:
Código
comparar_matrices <- function(P_teorica, P_empirica) {
n <- 100 # Suponemos 100 observaciones
esperada <- P_teorica * n
observada <- P_empirica * n
# Prueba chi cuadrado
test <- chisq.test(x = observada, p = esperada, simulate.p.value = TRUE)
return(list(test = test, P_empirica = P_empirica, P_teorica = P_teorica))
}Realizamos las comparaciones para cada k
Código
resultados_tests <- Map(comparar_matrices, lista_Pk, lista_Pk_empiricas)
resultados_tests[[1]]
[[1]]$test
Pearson's Chi-squared test with simulated p-value (based on 2000
replicates)
data: observada
X-squared = 7.011, df = NA, p-value = 0.01149
[[1]]$P_empirica
0 1
0 0.5172414 0.4827586
1 0.7000000 0.3000000
[[1]]$P_teorica
0 1
0 0.5893452 0.4106548
1 0.5954495 0.4045505
[[2]]
[[2]]$test
Pearson's Chi-squared test with simulated p-value (based on 2000
replicates)
data: observada
X-squared = 0.010662, df = NA, p-value = 1
[[2]]$P_empirica
0 1
0 0.6071429 0.3928571
1 0.6000000 0.4000000
[[2]]$P_teorica
0 1
0 0.6054697 0.3945303
1 0.5720690 0.4279310
[[3]]
[[3]]$test
Pearson's Chi-squared test with simulated p-value (based on 2000
replicates)
data: observada
X-squared = 0.18536, df = NA, p-value = 0.6742
[[3]]$P_empirica
0 1
0 0.6296296 0.3703704
1 0.6000000 0.4000000
[[3]]$P_teorica
0 1
0 0.5893452 0.4106548
1 0.5954495 0.4045505
[[4]]
[[4]]$test
Pearson's Chi-squared test with simulated p-value (based on 2000
replicates)
data: observada
X-squared = 0.04964, df = NA, p-value = 0.7836
[[4]]$P_empirica
0 1
0 0.6153846 0.3846154
1 0.6000000 0.4000000
[[4]]$P_teorica
0 1
0 0.5922921 0.4077079
1 0.5911765 0.4088235De Los resultados obtenidos podemos decir:
Para k = 1 , se observa una diferencia estadísticamente significativa (p-valor = 0.01099) entre las matrices empírica y teórica. Esto sugiere que el modelo teórico 1 puede no ajustarse perfectamente a los datos observados en el primer paso, indicando posiblemente la influencia de factores no modelados o variaciones aleatorias que no se capturan completamente en el modelo teórico inicial.
Para k = 2, 3, y 4, los p-valores son muy altos (1, 0.6752, y 0.7671 respectivamente), indicando que no hay diferencias estadísticamente significativas entre las matrices empíricas y teóricas en estos pasos. Esto implica que, a medida que avanzamos en más pasos, el modelo teórico parece capturar adecuadamente la dinámica del sistema según lo observado en los datos. Las pequeñas discrepancias observadas en \(k = 1\) podrían amortiguarse o distribuirse a través de múltiples transiciones, alineándose mejor las predicciones del modelo con las observaciones a medida que se consideran más pasos.
Estos resultados sugieren que, aunque puede haber algunas desviaciones iniciales entre las observaciones y las predicciones del modelo en el corto plazo (como en \(k = 1\)), el modelo de Markov \(P^{(1)}\) es efectivamente robusto en el análisis a más largo plazo, proporcionando una buena aproximación de las transiciones de estado a través de múltiples pasos.