Muestreo no probabilístico intencional

El investigador elige los casos que considera más representativos o informativos.
Se usa en estudios cualitativos o en investigaciones donde se busca un perfil específico.

Muestra

Participantes

De la poblacion total de 3.102 pacientes que ingresaron con traumatismo craneoencefálico(TCE) al Hospital Universitario Hernando Moncaleano Perdomo de Neiva, se tomo como muestra unicamente los pacientes con traumatismo moderado y severo, que participan en el Estudio “Secuelas cognitivas, comportamentales y emocionales con relación a la funcionalidad endocrina en el paciente con TCE ocasionado por accidentes de tránsito”, este estudio complementa la limitación respecto a la validación en esta condición clínica.

Un total de 104 participantes: - 34 con diagnóstico de TCE y 70 controles cognitivamente sanos, se tuvieron en cuenta los siguientes criterios clínicos:

Criterios de exclusión para los dos grupos TCE y controles:

  • Mujeres en estado de gestación.

  • Consumo activo de sustancias psicoactivas (SPA).

  • Historial neurológico y/o psiquiátrico.

  • Comorbilidades (diabetes, hipertensión arterial, dislipidemia) y/o alteraciones hormonales.

  • Historia clínica pre-mórbida de aprendizaje (se tendrá en cuenta la resistencia escolar).

  • Trastornos perceptivos visuales, auditivos o motores que imposibilitan el desarrollo de la evaluación.

Criterios de inclusión grupo clínico TCE - Hombres y mujeres entre los 18 a 50 años de edad.

  • Traumatismo craneoencefálico (TCE) secundario a accidente de tránsito, con puntuación en la escala de Glasgow en el rango moderado a severo.

  • Criterios de inclusión grupo control

  • Hombres y mujeres entre los 18 a 50 años de edad.

Cargar de Base de datos

# Leer el archivo datos.csv
datos <- read.csv("datos.csv", 
                  sep = ";",              
                  fileEncoding = "latin1", 
                  header = TRUE)          

Estimadores

Es importante dejar claro que, para este ejercicio únicamente haremos uso del grupo clínico.

TCE = 34 participantes, grupo clínico con TCE moderado y severo

Variables a trabajar:

1. Media Poblacional (μ):

edad <- datos$edad
media_estimada <- mean(edad)
media_estimada
## [1] 29.5
escolaridad <- datos$escolaridad
media_estimada <- mean(escolaridad)
media_estimada
## [1] 8.794118

La edad media de los participantes con traumatismo craneoencefálico y su promedio de años de escolaridad se muestran en los resultados anteriores.

2. Varianza Poblacional (σ2)

varianza_estimada <- var(edad)
desviacion_estandar <- sd(edad)

cat("Varianza de la edad:", varianza_estimada, "años²\n")
## Varianza de la edad: 87.04545 años²
cat("Desviación estándar:", desviacion_estandar, "años\n")
## Desviación estándar: 9.329815 años

La varianza de la edad muestra qué tan dispersos están los datos respecto a la media. La desviación estándar indica que, en promedio, las edades de los participantes se alejan de la media en aproximadamente la cantidad mostrada arriba.

3. Proporción Poblacional (p)

proporcion_por_Sexo <- table(datos$sexo) / nrow(datos)
proporcion_por_Sexo
## 
##  Femenino Masculino 
## 0.2058824 0.7941176

De acuerdo con los resultados, la mayoría de los participantes con TCE son de sexo masculino, mientras que una menor proporción corresponde al sexo femenino.

proporcion_por_Clasificación_TCE <- table(datos$clasificacion_tce) / nrow(datos)
proporcion_por_Clasificación_TCE
## 
##  Moderado    Severo 
## 0.5588235 0.4411765

De los participantes con TCE, una mayor proporción fueron clasificados como casos moderados, mientras que el resto correspondieron a casos severos, según la Escala de Coma de Glasgow (GCS).

proporcion_declive_moderado <- sum(datos$clasificacion_tce == "Moderado" & 
                                   datos$deterioro_clinico == "S?") /
                               sum(datos$clasificacion_tce == "Moderado")

proporcion_declive_moderado
## [1] 0

Algunos pacientes con TCE moderado presentaron deterioro clínico y pasaron a clasificarse como TCE severo, según la GCS.

4. Diferencia de Medias (μ1−μ2)

media_moderado <- mean(datos$edad[datos$clasificacion_tce == "Moderado"], na.rm = TRUE)
media_severo   <- mean(datos$edad[datos$clasificacion_tce == "Severo"], na.rm = TRUE)
diferencia_medias <- media_moderado - media_severo
media_moderado
## [1] 28.78947
media_severo
## [1] 30.4
diferencia_medias
## [1] -1.610526

Los resultados muestran las edades promedio de los pacientes con TCE moderado y severo. La diferencia de edad entre ambos grupos es pequeña, siendo los casos severos ligeramente mayores que los moderados.

media_moderado <- mean(datos$escolaridad
[datos$clasificacion_tce == "Moderado"], na.rm = TRUE)
media_severo   <- mean(datos$escolaridad
[datos$clasificacion_tce == "Severo"], na.rm = TRUE)
diferencia_medias <- media_moderado - media_severo
media_moderado
## [1] 9.263158
media_severo
## [1] 8.2
diferencia_medias
## [1] 1.063158

Los resultados muestran que los pacientes con TCE moderado tienen, en promedio, más años de escolaridad que aquellos con TCE severo. Esta diferencia podría sugerir una posible relación entre el nivel educativo y la gravedad del TCE.

5. Diferencia de Proporciones (p1−p2)

hombres_moderado <- sum(datos$clasificacion_tce == "Moderado" & datos$sexo == "Masculino", na.rm = TRUE)
hombres_severo   <- sum(datos$clasificacion_tce == "Severo"   & datos$sexo == "Masculino", na.rm = TRUE)
total_moderado <- sum(datos$clasificacion_tce == "Moderado", na.rm = TRUE)
total_severo   <- sum(datos$clasificacion_tce == "Severo", na.rm = TRUE)
p1 <- hombres_moderado / total_moderado
p2 <- hombres_severo   / total_severo
diferencia_p <- p1 - p2

# Mostrar proporciones
cat("Proporción de hombres en TCE moderado:", p1, "\n")
## Proporción de hombres en TCE moderado: 0.7368421
cat("Proporción de hombres en TCE severo:", p2, "\n")
## Proporción de hombres en TCE severo: 0.8666667
cat("Diferencia de proporciones (p1 - p2):", diferencia_p, "\n\n")
## Diferencia de proporciones (p1 - p2): -0.1298246
# Crear tabla de contingencia
tabla <- matrix(c(hombres_moderado, total_moderado - hombres_moderado,
                  hombres_severo, total_severo - hombres_severo),
                nrow = 2, byrow = TRUE,
                dimnames = list(c("Moderado", "Severo"), c("Masculino", "Femenino")))

print(tabla)
##          Masculino Femenino
## Moderado        14        5
## Severo          13        2
# Test Exacto de Fisher (más apropiado para muestras pequeñas)
fisher.test(tabla)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tabla
## p-value = 0.4263
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.03600613 3.30318204
## sample estimates:
## odds ratio 
##  0.4411833

En el análisis de la distribución por sexo, se observó que la proporción de hombres fue mayor entre los casos severos que en los moderados. Sin embargo, el test estadístico indica que esta diferencia no es significativa. Esto significa que, en esta muestra, la distribución por sexo es similar entre los pacientes con TCE moderado y severo.

Propiedades de estimadores

Variable: Estado_Civil

estado_civil <- as.factor(datos$estado_civil)
table(estado_civil)
## estado_civil
##     Casados    Solteros Unión_libre 
##           3          21          10

1. Sesgo

Un estimador insesgado de una proporción p es la proporción muestral (p̂ = x/n)

prop_civil <- prop.table(table(estado_civil))
cat("\n=== INSESGADEZ ===\n")
## 
## === INSESGADEZ ===
cat("Proporciones muestrales por categoría:\n")
## Proporciones muestrales por categoría:
print(prop_civil)
## estado_civil
##     Casados    Solteros Unión_libre 
##  0.08823529  0.61764706  0.29411765
cat("→ Las proporciones muestrales son estimadores insesgados de las proporciones poblacionales.\n")
## → Las proporciones muestrales son estimadores insesgados de las proporciones poblacionales.

Las proporciones calculadas de Estado_Civil no tienen sesgo, es decir, no sobreestiman ni subestiman las proporciones reales de la población. Los resultados de la muestra reflejan correctamente la realidad.

2. Eficiencia

En variables categóricas, la eficiencia se evalúa por la varianza del estimador de proporción: Var(p̂) = [p * (1 - p)] / n

n <- length(estado_civil)
varianza_prop <- prop_civil * (1 - prop_civil) / n

cat("\n=== EFICIENCIA ===\n")
## 
## === EFICIENCIA ===
cat("Varianza de los estimadores de proporción (menor = más eficiente):\n")
## Varianza de los estimadores de proporción (menor = más eficiente):
print(varianza_prop)
## estado_civil
##     Casados    Solteros Unión_libre 
## 0.002366171 0.006945858 0.006106249
cat("→ Las categorías con proporciones extremas tienen menor varianza.\n")
## → Las categorías con proporciones extremas tienen menor varianza.

Las categorías con proporciones más altas o más bajas tienen menor varianza, lo que significa que son estimaciones más precisas y estables.

3. Consistencia

La proporción muestral es consistente: al aumentar n, se aproxima a p real.

set.seed(123)
n_vals <- seq(10, n, by = 10)
prop_evol <- sapply(n_vals, function(k) {
  muestra <- sample(estado_civil, k, replace = TRUE)
  prop.table(table(muestra))["Solteros"]
})

plot(n_vals, prop_evol, type = "b", pch = 19, col = "darkblue",
     main = "Consistencia del estimador de proporción (Solteros)",
     xlab = "Tamaño de muestra", ylab = "Proporción estimada de Solteros")
abline(h = prop_civil["Solteros"], col = "red", lty = 2)

cat("\n=== CONSISTENCIA ===\n")
## 
## === CONSISTENCIA ===
cat("Observa el gráfico: la proporción estimada de Solteros se estabiliza conforme aumenta n.\n")
## Observa el gráfico: la proporción estimada de Solteros se estabiliza conforme aumenta n.

El gráfico muestra que cuando aumenta el tamaño de la muestra, la proporción estimada se acerca al valor real. Esto demuestra que el estimador es consistente.

4. Suficiencia

En una distribución binomial/multinomial, el conteo de éxitos (o frecuencias) es un estimador suficiente para p.

cat("\n=== SUFICIENCIA ===\n")
## 
## === SUFICIENCIA ===
cat("Para una variable categórica (multinomial), las frecuencias observadas son suficientes para estimar las proporciones poblacionales.\n")
## Para una variable categórica (multinomial), las frecuencias observadas son suficientes para estimar las proporciones poblacionales.
cat("→ La tabla de frecuencias contiene toda la información necesaria sobre p.\n")
## → La tabla de frecuencias contiene toda la información necesaria sobre p.

La tabla de frecuencias contiene toda la información necesaria para estimar las proporciones. No necesitamos datos adicionales para hacer estas estimaciones.

5. Robustez

Para categorías, la robustez se interpreta como estabilidad frente a errores de clasificación.

cat("\n=== Robustez ===\n")
## 
## === Robustez ===
cat("La proporción muestral es moderadamente robusta.\n")
## La proporción muestral es moderadamente robusta.
cat("Si existen pocos errores de clasificación en 'Estado_Civil', las proporciones cambian poco.\n")
## Si existen pocos errores de clasificación en 'Estado_Civil', las proporciones cambian poco.

El estimador de proporciones es robusto cuando pequeños errores en la clasificación no cambian mucho los resultados.

Simulación de cambiar aleatoriamente un 5% de respuestas

set.seed(123)
estado_mod <- estado_civil
indices <- sample(1:n, size = 0.05*n)
estado_mod[indices] <- sample(levels(estado_civil), length(indices), replace = TRUE)

prop_original <- prop.table(table(estado_civil))
prop_modificada <- prop.table(table(estado_mod))

cat("\nProporciones originales:\n")
## 
## Proporciones originales:
print(prop_original)
## estado_civil
##     Casados    Solteros Unión_libre 
##  0.08823529  0.61764706  0.29411765
cat("Proporciones con 5% de errores:\n")
## Proporciones con 5% de errores:
print(prop_modificada)
## estado_mod
##     Casados    Solteros Unión_libre 
##  0.08823529  0.61764706  0.29411765
cat("→ Cambios pequeños indican robustez del estimador.\n")
## → Cambios pequeños indican robustez del estimador.

Al introducir 5% de errores aleatorios, se puede observar si el estimador es robusto. Si las proporciones cambian poco, el estimador es robusto; si cambian mucho, no lo es. En este caso, los cambios observados permiten evaluar la estabilidad del estimador frente a errores de clasificación.

Intervalos de Confianza

Un intervalo de confianza (IC) es un rango de valores construido a partir de una muestra que se utiliza para estimar un parámetro poblacional desconocido. Incorpora la variabilidad inherente al muestreo.

La estimación puntual (e.g., la media muestral) no informa cuánta incertidumbre existe. El IC del 95% incorpora esa incertidumbre: si repitiéramos el estudio muchas veces, aproximadamente el 95% de los intervalos construidos contendrían el verdadero parámetro poblacional.

Variable 1: Escolaridad

Para esta sección, analizaremos la escolaridad (años de educación formal) de los pacientes con TCE, construyendo un intervalo de confianza del 95% para la media poblacional.

Estadísticas Descriptivas

# Extraer la variable escolaridad
escolaridad_vec <- datos$escolaridad

# Calcular estadísticas descriptivas
n_escol <- length(escolaridad_vec)
media_escol <- mean(escolaridad_vec, na.rm = TRUE)
sd_escol <- sd(escolaridad_vec, na.rm = TRUE)
min_escol <- min(escolaridad_vec, na.rm = TRUE)
max_escol <- max(escolaridad_vec, na.rm = TRUE)

cat("=== ESTADÍSTICAS DESCRIPTIVAS: Escolaridad ===\n")
## === ESTADÍSTICAS DESCRIPTIVAS: Escolaridad ===
cat("Tamaño de muestra:", n_escol, "\n")
## Tamaño de muestra: 34
cat("Media:", round(media_escol, 2), "años\n")
## Media: 8.79 años
cat("Desviación estándar:", round(sd_escol, 2), "años\n")
## Desviación estándar: 4.15 años
cat("Mínimo:", min_escol, "años\n")
## Mínimo: 2 años
cat("Máximo:", max_escol, "años\n")
## Máximo: 16 años

Los resultados muestran el promedio de años de escolaridad de los pacientes con TCE, así como la variabilidad de estos datos.

Intervalo de Confianza del 95% para la Media

# Nivel de confianza del 95%
nivel_confianza <- 0.95
alpha <- 1 - nivel_confianza

# Error estándar de la media
error_estandar <- sd_escol / sqrt(n_escol)

# Valor crítico de la distribución t-Student
valor_critico <- qt(1 - alpha/2, df = n_escol - 1)

# Margen de error
margen_error <- valor_critico * error_estandar

# Límites del intervalo de confianza
limite_inferior <- media_escol - margen_error
limite_superior <- media_escol + margen_error

cat("\n=== INTERVALO DE CONFIANZA DEL 95% PARA LA MEDIA ===\n")
## 
## === INTERVALO DE CONFIANZA DEL 95% PARA LA MEDIA ===
cat("Media muestral:", round(media_escol, 2), "años\n")
## Media muestral: 8.79 años
cat("Error estándar:", round(error_estandar, 2), "\n")
## Error estándar: 0.71
cat("Valor crítico (t):", round(valor_critico, 3), "\n")
## Valor crítico (t): 2.035
cat("Margen de error:", round(margen_error, 2), "años\n")
## Margen de error: 1.45 años
cat("IC 95%: [", round(limite_inferior, 2), ",", round(limite_superior, 2), "] años\n")
## IC 95%: [ 7.35 , 10.24 ] años

El IC del 95% para la escolaridad indica el rango plausible de la media poblacional de años de educación formal en pacientes con TCE moderado y severo. Este dato es clínicamente relevante porque un mayor nivel educativo se asocia con mayor reserva cognitiva, lo que puede modular la recuperación neuropsicológica tras la lesión.

Simulación de Múltiples Intervalos de Confianza

Para comprender mejor la interpretación del IC, realizaremos una simulación que muestra cómo diferentes muestras generan diferentes intervalos de confianza.

set.seed(456)

# Parámetros de la simulación
n_muestra <- n_escol        # tamaño de cada muestra (34 pacientes)
mu_estimado <- media_escol  # usamos la media muestral como estimación de la media poblacional
sigma_estimado <- sd_escol  # usamos la desviación estándar muestral

# Número de simulaciones a realizar
Nsim <- 15  # número de intervalos a simular

cat("Número de simulaciones a realizar:", Nsim, "\n")
## Número de simulaciones a realizar: 15
# Crear dataframe para almacenar resultados
resultados_sim <- data.frame(
  iteracion = 1:Nsim,
  limite_inf = NA,
  limite_sup = NA,
  contiene_media = NA
)

# Simulación de múltiples muestras e intervalos
for (i in 1:Nsim) {
  # Generar una muestra aleatoria (asumiendo distribución normal)
  muestra_sim <- rnorm(n_muestra, mean = mu_estimado, sd = sigma_estimado)
  
  # Calcular media y desviación estándar de la muestra
  media_sim <- mean(muestra_sim)
  sd_sim <- sd(muestra_sim)
  se_sim <- sd_sim / sqrt(n_muestra)
  
  # Valor crítico t
  t_crit <- qt(0.975, df = n_muestra - 1)
  
  # Calcular límites del IC
  L <- media_sim - t_crit * se_sim
  U <- media_sim + t_crit * se_sim
  
  # Guardar resultados
  resultados_sim$limite_inf[i] <- L
  resultados_sim$limite_sup[i] <- U
  resultados_sim$contiene_media[i] <- (L <= mu_estimado & U >= mu_estimado)
}

# Calcular porcentaje de intervalos que contienen la media
porcentaje_cobertura <- sum(resultados_sim$contiene_media) / Nsim * 100

cat("\n=== RESULTADOS DE LA SIMULACIÓN ===\n")
## 
## === RESULTADOS DE LA SIMULACIÓN ===
cat("Número de intervalos simulados:", Nsim, "\n")
## Número de intervalos simulados: 15
cat("Intervalos que contienen la media:", sum(resultados_sim$contiene_media), "\n")
## Intervalos que contienen la media: 14
cat("Porcentaje de cobertura:", round(porcentaje_cobertura, 1), "%\n")
## Porcentaje de cobertura: 93.3 %

Visualización de la Simulación

# Crear gráfico de los intervalos de confianza
plot(NA, 
     xlim = c(min(resultados_sim$limite_inf), max(resultados_sim$limite_sup)),
     ylim = c(1, Nsim),
     xlab = "Años de Escolaridad", 
     ylab = "Número de Intervalo",
     main = paste0(Nsim, " Intervalos de Confianza del 95% para Escolaridad"))

# Dibujar cada intervalo
for (i in 1:Nsim) {
  color_linea <- ifelse(resultados_sim$contiene_media[i], "blue", "red")
  segments(resultados_sim$limite_inf[i], i, 
           resultados_sim$limite_sup[i], i, 
           col = color_linea, lwd = 2)
}

# Línea vertical en la media estimada
abline(v = mu_estimado, col = "darkgreen", lwd = 3, lty = 2)

# Leyenda
legend("topright", 
       legend = c("Contiene la media", "No contiene la media", "Media estimada"),
       col = c("blue", "red", "darkgreen"), 
       lwd = c(2, 2, 3),
       lty = c(1, 1, 2),
       cex = 0.8)

El gráfico muestra múltiples intervalos de confianza del 95%, cada uno calculado a partir de muestras simuladas:

  • Líneas azules: intervalos que SÍ contienen la media estimada de la población
  • Líneas rojas: intervalos que NO contienen la media estimada
  • Línea verde vertical: la media estimada de años de escolaridad (valor de referencia)

Este gráfico ilustra que, aunque construimos intervalos con 95% de confianza, algunos intervalos (aproximadamente el 5%) no contendrán el verdadero valor del parámetro poblacional. Sin embargo, en la práctica, solo disponemos de un intervalo (el calculado con nuestra muestra real) y no sabemos si es uno de los que contiene o no el valor verdadero.

Comparación del IC por Clasificación de TCE

Finalmente, comparemos los intervalos de confianza de escolaridad entre pacientes con TCE moderado y severo.

# Separar datos por clasificación
escol_moderado <- datos$escolaridad[datos$clasificacion_tce == "Moderado"]
escol_severo <- datos$escolaridad[datos$clasificacion_tce == "Severo"]

# Función para calcular IC del 95%
calcular_ic <- function(x, nivel = 0.95) {
  n <- length(x)
  media <- mean(x, na.rm = TRUE)
  sd_x <- sd(x, na.rm = TRUE)
  se <- sd_x / sqrt(n)
  t_crit <- qt((1 + nivel)/2, df = n - 1)
  margen <- t_crit * se
  return(c(media = media, 
           limite_inf = media - margen, 
           limite_sup = media + margen,
           n = n))
}

# Calcular IC para ambos grupos
ic_moderado <- calcular_ic(escol_moderado)
ic_severo <- calcular_ic(escol_severo)

cat("\n=== INTERVALOS DE CONFIANZA POR CLASIFICACIÓN TCE ===\n\n")
## 
## === INTERVALOS DE CONFIANZA POR CLASIFICACIÓN TCE ===
cat("TCE MODERADO:\n")
## TCE MODERADO:
cat("  n =", ic_moderado["n"], "\n")
##   n = 19
cat("  Media =", round(ic_moderado["media"], 2), "años\n")
##   Media = 9.26 años
cat("  IC 95% = [", round(ic_moderado["limite_inf"], 2), ",", 
    round(ic_moderado["limite_sup"], 2), "] años\n\n")
##   IC 95% = [ 7.3 , 11.23 ] años
cat("TCE SEVERO:\n")
## TCE SEVERO:
cat("  n =", ic_severo["n"], "\n")
##   n = 15
cat("  Media =", round(ic_severo["media"], 2), "años\n")
##   Media = 8.2 años
cat("  IC 95% = [", round(ic_severo["limite_inf"], 2), ",", 
    round(ic_severo["limite_sup"], 2), "] años\n")
##   IC 95% = [ 5.82 , 10.58 ] años

Los resultados permiten comparar el nivel educativo entre pacientes con TCE moderado y severo. Las diferencias en escolaridad pueden influir en el pronóstico y rehabilitación cognitiva.

Visualización Comparativa

# Crear gráfico comparativo
par(mar = c(5, 6, 4, 2))

plot(1:2, c(ic_moderado["media"], ic_severo["media"]), 
     xlim = c(0.5, 2.5), 
     ylim = c(0, max(c(ic_moderado["limite_sup"], ic_severo["limite_sup"])) * 1.1),
     xlab = "", ylab = "Años de Escolaridad", 
     main = "Comparación de IC 95% de Escolaridad por Clasificación de TCE",
     pch = 19, cex = 2, col = c("steelblue", "coral"),
     xaxt = "n")

# Añadir intervalos de confianza
arrows(1, ic_moderado["limite_inf"], 
       1, ic_moderado["limite_sup"], 
       angle = 90, code = 3, length = 0.1, lwd = 2, col = "steelblue")
arrows(2, ic_severo["limite_inf"], 
       2, ic_severo["limite_sup"], 
       angle = 90, code = 3, length = 0.1, lwd = 2, col = "coral")

# Etiquetas del eje x
axis(1, at = 1:2, labels = c("TCE Moderado", "TCE Severo"))

# Añadir valores de las medias
text(1, ic_moderado["media"], 
     labels = round(ic_moderado["media"], 1), 
     pos = 4, cex = 0.9)
text(2, ic_severo["media"], 
     labels = round(ic_severo["media"], 1), 
     pos = 4, cex = 0.9)

# Línea de referencia en cero
abline(h = 0, lty = 3, col = "gray")

El gráfico muestra las medias y sus intervalos de confianza del 95% para ambos grupos. Los puntos representan las medias muestrales y las barras verticales representan los intervalos de confianza. Si los intervalos no se solapan, esto sugiere una diferencia significativa en los años de escolaridad entre ambos grupos.

Variable 2: ACER Total (Puntuación Cognitiva)

Ahora analizaremos el ACER Total (Addenbrooke’s Cognitive Examination-Revised), que es una prueba cognitiva estandarizada con puntuación de 0 a 100 puntos. Esta variable es fundamental para evaluar el estado cognitivo de los pacientes con TCE, siendo más sensible que el MMSE para detectar deterioro cognitivo leve.

Estadísticas Descriptivas del ACER

# Extraer la variable ACER total
acer_vec <- datos$hacer_total_100

# Calcular estadísticas descriptivas
n_acer <- length(acer_vec)
media_acer <- mean(acer_vec, na.rm = TRUE)
sd_acer <- sd(acer_vec, na.rm = TRUE)
min_acer <- min(acer_vec, na.rm = TRUE)
max_acer <- max(acer_vec, na.rm = TRUE)

cat("=== ESTADÍSTICAS DESCRIPTIVAS: ACER Total (sobre 100 puntos) ===\n")
## === ESTADÍSTICAS DESCRIPTIVAS: ACER Total (sobre 100 puntos) ===
cat("Tamaño de muestra:", n_acer, "\n")
## Tamaño de muestra: 34
cat("Media:", round(media_acer, 2), "puntos\n")
## Media: 62.79 puntos
cat("Desviación estándar:", round(sd_acer, 2), "puntos\n")
## Desviación estándar: 19.67 puntos
cat("Mínimo:", min_acer, "puntos\n")
## Mínimo: 15 puntos
cat("Máximo:", max_acer, "puntos\n")
## Máximo: 92 puntos

El ACER es una herramienta de evaluación cognitiva ampliamente utilizada en neuropsicología. Puntuaciones más bajas indican mayor deterioro cognitivo. Un punto de corte típico es 88 puntos para detectar deterioro cognitivo.

Intervalo de Confianza del 95% para el ACER

# Nivel de confianza del 95%
nivel_confianza <- 0.95
alpha <- 1 - nivel_confianza

# Error estándar de la media
error_estandar_acer <- sd_acer / sqrt(n_acer)

# Valor crítico de la distribución normal estándar (Z)
valor_critico_acer <- qnorm(1 - alpha/2)

# Margen de error
margen_error_acer <- valor_critico_acer * error_estandar_acer

# Límites del intervalo de confianza
limite_inferior_acer <- media_acer - margen_error_acer
limite_superior_acer <- media_acer + margen_error_acer

cat("\n=== INTERVALO DE CONFIANZA DEL 95% PARA ACER ===\n")
## 
## === INTERVALO DE CONFIANZA DEL 95% PARA ACER ===
cat("Media muestral:", round(media_acer, 2), "puntos\n")
## Media muestral: 62.79 puntos
cat("Error estándar:", round(error_estandar_acer, 2), "\n")
## Error estándar: 3.37
cat("Valor crítico (Z):", round(valor_critico_acer, 4), "\n")
## Valor crítico (Z): 1.96
cat("Margen de error:", round(margen_error_acer, 2), "puntos\n")
## Margen de error: 6.61 puntos
cat("IC 95%: [", round(limite_inferior_acer, 2), ",", round(limite_superior_acer, 2), "] puntos\n")
## IC 95%: [ 56.18 , 69.4 ] puntos

El IC del 95% para el ACER delimita el rango plausible del rendimiento cognitivo medio en la población TCE. Comparar este intervalo con el punto de corte clínico de 88 puntos permite determinar si el grupo, en conjunto, se encuentra en zona de deterioro cognitivo probable.

Comparación del ACER por Clasificación de TCE

Comparemos el rendimiento cognitivo (ACER) entre pacientes con TCE moderado y severo.

# Separar datos por clasificación
acer_moderado <- datos$hacer_total_100[datos$clasificacion_tce == "Moderado"]
acer_severo <- datos$hacer_total_100[datos$clasificacion_tce == "Severo"]

# Función para calcular IC del 95% usando distribución normal (Z)
calcular_ic_z <- function(x, nivel = 0.95) {
  n <- length(x)
  media <- mean(x, na.rm = TRUE)
  sd_x <- sd(x, na.rm = TRUE)
  se <- sd_x / sqrt(n)
  z_crit <- qnorm((1 + nivel)/2)
  margen <- z_crit * se
  return(c(media = media, 
           limite_inf = media - margen, 
           limite_sup = media + margen,
           n = n))
}

# Calcular IC para ambos grupos
ic_acer_moderado <- calcular_ic_z(acer_moderado)
ic_acer_severo <- calcular_ic_z(acer_severo)

cat("\n=== INTERVALOS DE CONFIANZA ACER POR CLASIFICACIÓN TCE ===\n\n")
## 
## === INTERVALOS DE CONFIANZA ACER POR CLASIFICACIÓN TCE ===
cat("TCE MODERADO:\n")
## TCE MODERADO:
cat("  n =", ic_acer_moderado["n"], "\n")
##   n = 19
cat("  Media =", round(ic_acer_moderado["media"], 2), "puntos\n")
##   Media = 67.32 puntos
cat("  IC 95% = [", round(ic_acer_moderado["limite_inf"], 2), ",", 
    round(ic_acer_moderado["limite_sup"], 2), "] puntos\n\n")
##   IC 95% = [ 59.73 , 74.9 ] puntos
cat("TCE SEVERO:\n")
## TCE SEVERO:
cat("  n =", ic_acer_severo["n"], "\n")
##   n = 15
cat("  Media =", round(ic_acer_severo["media"], 2), "puntos\n")
##   Media = 57.07 puntos
cat("  IC 95% = [", round(ic_acer_severo["limite_inf"], 2), ",", 
    round(ic_acer_severo["limite_sup"], 2), "] puntos\n")
##   IC 95% = [ 45.95 , 68.19 ] puntos

Los resultados permiten comparar el rendimiento cognitivo entre ambos grupos. Se espera que los pacientes con TCE severo presenten puntuaciones más bajas en el ACER, reflejando mayor deterioro cognitivo.

Visualización Comparativa del ACER

# Crear gráfico comparativo
par(mar = c(5, 6, 4, 2))

plot(1:2, c(ic_acer_moderado["media"], ic_acer_severo["media"]), 
     xlim = c(0.5, 2.5), 
     ylim = c(0, 100),  # Escala del ACER es 0-100
     xlab = "", ylab = "ACER Total (puntos)", 
     main = "Comparación de IC 95% de ACER por Clasificación de TCE",
     pch = 19, cex = 2, col = c("steelblue", "coral"),
     xaxt = "n")

# Añadir intervalos de confianza
arrows(1, ic_acer_moderado["limite_inf"], 
       1, ic_acer_moderado["limite_sup"], 
       angle = 90, code = 3, length = 0.1, lwd = 2, col = "steelblue")
arrows(2, ic_acer_severo["limite_inf"], 
       2, ic_acer_severo["limite_sup"], 
       angle = 90, code = 3, length = 0.1, lwd = 2, col = "coral")

# Etiquetas del eje x
axis(1, at = 1:2, labels = c("TCE Moderado", "TCE Severo"))

# Añadir valores de las medias
text(1, ic_acer_moderado["media"], 
     labels = round(ic_acer_moderado["media"], 1), 
     pos = 4, cex = 0.9)
text(2, ic_acer_severo["media"], 
     labels = round(ic_acer_severo["media"], 1), 
     pos = 4, cex = 0.9)

# Línea de referencia en 88 (punto de corte común para deterioro cognitivo)
abline(h = 88, lty = 2, col = "red", lwd = 1.5)
text(2.3, 88, "Punto de corte (88)", pos = 3, cex = 0.8, col = "red")

# Línea de referencia en cero
abline(h = 0, lty = 3, col = "gray")

El gráfico compara las puntuaciones del ACER entre ambos grupos. La línea roja marca el punto de corte típico (88 puntos) que se utiliza para identificar deterioro cognitivo. Puntuaciones por debajo de 88 sugieren deterioro. Las diferencias entre grupos son relevantes para entender el impacto de la severidad del TCE en la función cognitiva.

Prueba de Hipótesis: Correlación entre Nivel Escolar y ACE-R

Planteamiento de Hipótesis

Hipótesis Nula (H₀): No existe una relación estadísticamente significativa entre el nivel escolar y los resultados obtenidos en el Addenbrooke’s Cognitive Examination Revised (ACE-R) en pacientes con Traumatismo Craneoencefálico (TCE).

  • En términos estadísticos: ρ = 0

Hipótesis Alternativa (H₁): Existe una relación positiva estadísticamente significativa entre el nivel escolar y los resultados obtenidos en el ACE-R en pacientes con TCE. Es decir, los pacientes con mayor nivel escolar obtendrán, en promedio, mayores puntuaciones en el ACE-R.

  • En términos estadísticos: ρ > 0

Análisis Exploratorio

Primero, examinaremos visualmente la relación entre ambas variables mediante un diagrama de dispersión.

# Extraer las variables
escolaridad_test <- datos$escolaridad
acer_test <- datos$hacer_total_100

# Crear diagrama de dispersión
plot(escolaridad_test, acer_test,
     xlab = "Años de Escolaridad",
     ylab = "Puntuación ACE-R (sobre 100)",
     main = "Relación entre Escolaridad y ACE-R en Pacientes con TCE",
     pch = 19, col = "steelblue", cex = 1.5)

# Añadir línea de regresión
abline(lm(acer_test ~ escolaridad_test), col = "red", lwd = 2)

# Añadir línea de referencia del punto de corte ACE-R
abline(h = 88, lty = 2, col = "gray")
text(max(escolaridad_test) * 0.9, 88, "Punto de corte (88)", pos = 3, cex = 0.8)

# Calcular correlación preliminar
cor_preliminar <- cor(escolaridad_test, acer_test, use = "complete.obs")
text(min(escolaridad_test) + 2, max(acer_test) - 5, 
     paste0("r = ", round(cor_preliminar, 3)), 
     cex = 1.2, col = "darkred")

El diagrama de dispersión permite visualizar la tendencia de la relación entre ambas variables. La línea roja representa la tendencia lineal.

Verificación de Supuestos

Dado que el diagrama de dispersión muestra datos considerablemente dispersos con posibles valores atípicos, utilizaremos el test de correlación de Spearman en lugar de Pearson. Spearman es más robusto ante:

  • Datos dispersos y outliers
  • Violaciones de normalidad
  • Relaciones monótonas no necesariamente lineales
# Test de normalidad Shapiro-Wilk
shapiro_escolaridad <- shapiro.test(escolaridad_test)
shapiro_acer <- shapiro.test(acer_test)

cat("=== PRUEBAS DE NORMALIDAD ===\n\n")
## === PRUEBAS DE NORMALIDAD ===
cat("Escolaridad (Shapiro-Wilk):\n")
## Escolaridad (Shapiro-Wilk):
cat("  W =", round(shapiro_escolaridad$statistic, 4), "\n")
##   W = 0.9276
cat("  p-valor =", round(shapiro_escolaridad$p.value, 4), "\n")
##   p-valor = 0.0266
cat("  Interpretación:", ifelse(shapiro_escolaridad$p.value > 0.05, 
                                  "Los datos parecen seguir una distribución normal",
                                  "Los datos NO siguen una distribución normal"), "\n\n")
##   Interpretación: Los datos NO siguen una distribución normal
cat("ACE-R (Shapiro-Wilk):\n")
## ACE-R (Shapiro-Wilk):
cat("  W =", round(shapiro_acer$statistic, 4), "\n")
##   W = 0.9537
cat("  p-valor =", round(shapiro_acer$p.value, 4), "\n")
##   p-valor = 0.1582
cat("  Interpretación:", ifelse(shapiro_acer$p.value > 0.05, 
                                  "Los datos parecen seguir una distribución normal",
                                  "Los datos NO siguen una distribución normal"), "\n\n")
##   Interpretación: Los datos parecen seguir una distribución normal
cat("=== DECISIÓN: Usar correlación de Spearman ===\n")
## === DECISIÓN: Usar correlación de Spearman ===
cat("Debido a la dispersión de los datos observada en el gráfico y la presencia\n")
## Debido a la dispersión de los datos observada en el gráfico y la presencia
cat("de posibles valores atípicos, se utilizará el coeficiente de correlación\n")
## de posibles valores atípicos, se utilizará el coeficiente de correlación
cat("de Spearman (rho), que evalúa la correlación de rangos en lugar de valores\n")
## de Spearman (rho), que evalúa la correlación de rangos en lugar de valores
cat("directos, siendo más robusto ante estas condiciones.\n")
## directos, siendo más robusto ante estas condiciones.

Verificación visual de normalidad (QQ-plots)

par(mfrow = c(1, 2))

qqnorm(escolaridad_test,
       main = "QQ-Plot: Escolaridad",
       xlab = "Cuantiles teóricos", ylab = "Cuantiles observados",
       pch = 19, col = "steelblue", cex = 1.2)
qqline(escolaridad_test, col = "red", lwd = 2)

qqnorm(acer_test,
       main = "QQ-Plot: ACE-R",
       xlab = "Cuantiles teóricos", ylab = "Cuantiles observados",
       pch = 19, col = "steelblue", cex = 1.2)
qqline(acer_test, col = "red", lwd = 2)

par(mfrow = c(1, 1))

Los QQ-plots complementan el test de Shapiro-Wilk con una evaluación visual: los puntos deben alinearse sobre la línea roja si los datos siguen una distribución normal. Las desviaciones en los extremos —especialmente en la Escolaridad— confirman la presencia de asimetría o colas pesadas, lo que refuerza la elección de Spearman sobre Pearson.

Prueba de Correlación de Spearman (Una Cola)

Realizaremos la prueba de correlación de Spearman con una alternativa unilateral (cola derecha), ya que nuestra hipótesis alternativa establece que ρ > 0. Spearman evalúa la correlación basándose en los rangos de los datos, lo que lo hace más robusto ante datos dispersos y valores atípicos.

# Test de correlación de Spearman (una cola, alternativa "greater")
test_correlacion <- cor.test(escolaridad_test, acer_test, 
                              method = "spearman",
                              alternative = "greater",
                              exact = FALSE)

cat("\n=== PRUEBA DE CORRELACIÓN DE SPEARMAN (UNA COLA) ===\n\n")
## 
## === PRUEBA DE CORRELACIÓN DE SPEARMAN (UNA COLA) ===
cat("Coeficiente de correlación de rangos (rho):", round(test_correlacion$estimate, 4), "\n")
## Coeficiente de correlación de rangos (rho): 0.6451
cat("Estadístico S:", round(test_correlacion$statistic, 4), "\n")
## Estadístico S: 2322.863
cat("p-valor (una cola):", round(test_correlacion$p.value, 4), "\n\n")
## p-valor (una cola): 0
cat("Nota: Spearman evalúa la correlación entre los RANGOS de las variables,\n")
## Nota: Spearman evalúa la correlación entre los RANGOS de las variables,
cat("no los valores directos, lo que lo hace robusto ante datos dispersos.\n")
## no los valores directos, lo que lo hace robusto ante datos dispersos.
# Nivel de significancia
alpha <- 0.05

cat("=== DECISIÓN ESTADÍSTICA (α = 0.05) ===\n\n")
## === DECISIÓN ESTADÍSTICA (α = 0.05) ===
if (test_correlacion$p.value < alpha) {
  cat("✓ RECHAZAMOS H₀\n\n")
  cat("Conclusión: Existe evidencia estadística suficiente para afirmar que hay una\n")
  cat("relación monótona positiva significativa entre el nivel escolar y las puntuaciones\n")
  cat("del ACE-R en pacientes con TCE (rho =", round(test_correlacion$estimate, 3), 
      ", p =", round(test_correlacion$p.value, 4), ").\n\n")
  cat("Interpretación: A pesar de la dispersión en los datos, existe una tendencia\n")
  cat("significativa donde pacientes con mayor escolaridad tienden a obtener mejores\n")
  cat("puntuaciones en el ACE-R. Esto sugiere que la educación puede tener un efecto\n")
  cat("protector sobre el funcionamiento cognitivo en pacientes con TCE, aunque la\n")
  cat("variabilidad individual es considerable.\n")
} else {
  cat("✗ NO RECHAZAMOS H₀\n\n")
  cat("Conclusión: No existe evidencia estadística suficiente para afirmar que hay una\n")
  cat("relación monótona positiva significativa entre el nivel escolar y las puntuaciones\n")
  cat("del ACE-R en pacientes con TCE (rho =", round(test_correlacion$estimate, 3), 
      ", p =", round(test_correlacion$p.value, 4), ").\n\n")
  cat("Interpretación: Aunque pudiera existir cierta tendencia, la alta dispersión\n")
  cat("en los datos y la variabilidad individual no permiten concluir que el nivel\n")
  cat("educativo está significativamente relacionado con el rendimiento cognitivo en el ACE-R.\n")
}
## ✓ RECHAZAMOS H₀
## 
## Conclusión: Existe evidencia estadística suficiente para afirmar que hay una
## relación monótona positiva significativa entre el nivel escolar y las puntuaciones
## del ACE-R en pacientes con TCE (rho = 0.645 , p = 0 ).
## 
## Interpretación: A pesar de la dispersión en los datos, existe una tendencia
## significativa donde pacientes con mayor escolaridad tienden a obtener mejores
## puntuaciones en el ACE-R. Esto sugiere que la educación puede tener un efecto
## protector sobre el funcionamiento cognitivo en pacientes con TCE, aunque la
## variabilidad individual es considerable.

Interpretación del Coeficiente de Correlación

cat("\n=== MAGNITUD DE LA CORRELACIÓN ===\n\n")
## 
## === MAGNITUD DE LA CORRELACIÓN ===
r_valor <- abs(test_correlacion$estimate)

if (r_valor < 0.3) {
  magnitud <- "débil"
} else if (r_valor < 0.5) {
  magnitud <- "moderada"
} else if (r_valor < 0.7) {
  magnitud <- "fuerte"
} else {
  magnitud <- "muy fuerte"
}

cat("El coeficiente de correlación de Spearman (rho) =", round(test_correlacion$estimate, 3), 
    "indica una relación monótona", magnitud, "entre las variables.\n\n")
## El coeficiente de correlación de Spearman (rho) = 0.645 indica una relación monótona fuerte entre las variables.
cat("Nota importante: A diferencia del coeficiente de Pearson, el de Spearman evalúa\n")
## Nota importante: A diferencia del coeficiente de Pearson, el de Spearman evalúa
cat("la fuerza de la relación monótona (creciente o decreciente) entre los rangos\n")
## la fuerza de la relación monótona (creciente o decreciente) entre los rangos
cat("de las variables, sin asumir una relación lineal perfecta. Esto es especialmente\n")
## de las variables, sin asumir una relación lineal perfecta. Esto es especialmente
cat("útil cuando los datos presentan dispersión y posibles valores atípicos.\n")
## útil cuando los datos presentan dispersión y posibles valores atípicos.

Resumen de Resultados

cat("\n", rep("=", 70), "\n", sep = "")
## 
## ======================================================================
cat("                    RESUMEN DE LA PRUEBA DE HIPÓTESIS\n")
##                     RESUMEN DE LA PRUEBA DE HIPÓTESIS
cat(rep("=", 70), "\n\n", sep = "")
## ======================================================================
cat("Variables analizadas:\n")
## Variables analizadas:
cat("  • Variable independiente: Escolaridad (años de educación)\n")
##   • Variable independiente: Escolaridad (años de educación)
cat("  • Variable dependiente: ACE-R Total (puntuación cognitiva)\n\n")
##   • Variable dependiente: ACE-R Total (puntuación cognitiva)
cat("Tamaño de muestra: n =", nrow(na.omit(cbind(escolaridad_test, acer_test))), "\n\n")
## Tamaño de muestra: n = 34
cat("Prueba estadística: Correlación de Spearman (una cola)\n")
## Prueba estadística: Correlación de Spearman (una cola)
cat("Nivel de significancia: α = 0.05\n\n")
## Nivel de significancia: α = 0.05
cat("Resultados:\n")
## Resultados:
cat("  • rho (Spearman) =", round(test_correlacion$estimate, 4), "\n")
##   • rho (Spearman) = 0.6451
cat("  • p-valor =", round(test_correlacion$p.value, 4), "\n")
##   • p-valor = 0
cat("  • Decisión:", ifelse(test_correlacion$p.value < 0.05, 
                             "Se rechaza H₀", "No se rechaza H₀"), "\n\n")
##   • Decisión: Se rechaza H₀
cat(rep("=", 70), "\n", sep = "")
## ======================================================================

Los resultados de la prueba de Spearman permiten concluir si la escolaridad actúa como indicador de reserva cognitiva en esta población: una mayor educación formal suele asociarse con mayor densidad sináptica y plasticidad cerebral, factores que pueden atenuar el impacto funcional del TCE. La variabilidad individual observada en los datos sugiere, no obstante, que la gravedad de la lesión y otros factores clínicos (edad, localización, tiempo de evolución) también modulan el desempeño cognitivo post-TCE, lo que apunta a la necesidad de un abordaje rehabilitador individualizado.


Modelos de Regresión

1. Planteamiento del Problema e Hipótesis

En esta sección se busca identificar si algunas variables del estudio permiten explicar o predecir el rendimiento cognitivo de los pacientes con TCE, medido a través del ACE-R. Se trabaja con cuatro variables que, desde el punto de vista clínico, podrían tener relación con el desempeño cognitivo.

Variable respuesta (Y): hacer_total_100 — Puntuación total del ACE-R (0–100 puntos).

Variables predictoras:

Variable Descripción Escala
escolaridad Años de educación formal Continua
edad Edad del paciente (años) Continua
dias_hospitalizacion Días de estancia hospitalaria Continua
glasgow_final Puntuación Glasgow al egreso (3–15) Continua

Hipótesis:

  • H₀: Ninguna de las variables consideradas tiene una relación lineal significativa con el ACE-R. β₁ = β₂ = β₃ = β₄ = 0
  • H₁: Al menos una de las variables tiene una relación lineal significativa con el ACE-R.

2. Análisis de Correlación

Librerías adicionales

library(lmtest)
## Warning: package 'lmtest' was built under R version 4.5.2
## Cargando paquete requerido: zoo
## Warning: package 'zoo' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
if (!require("car", quietly = TRUE)) install.packages("car")
## Warning: package 'car' was built under R version 4.5.2
## Warning: package 'carData' was built under R version 4.5.2
library(car)

Estadísticas descriptivas de las variables del modelo

vars_modelo <- datos[, c("hacer_total_100", "escolaridad", "edad",
                          "dias_hospitalizacion", "glasgow_final")]
vars_modelo <- na.omit(vars_modelo)

cat(paste0(
  "=== ESTADÍSTICAS DESCRIPTIVAS ===\n\n",
  paste(capture.output(summary(vars_modelo)), collapse = "\n"),
  "\n\nN observaciones válidas: ", nrow(vars_modelo)
))
## === ESTADÍSTICAS DESCRIPTIVAS ===
## 
##  hacer_total_100  escolaridad          edad       dias_hospitalizacion
##  Min.   :15.00   Min.   : 2.000   Min.   :18.00   Min.   : 3.00       
##  1st Qu.:50.00   1st Qu.: 5.000   1st Qu.:21.25   1st Qu.: 9.00       
##  Median :64.50   Median :10.000   Median :27.50   Median :18.00       
##  Mean   :62.79   Mean   : 8.794   Mean   :29.50   Mean   :20.26       
##  3rd Qu.:77.75   3rd Qu.:12.000   3rd Qu.:36.75   3rd Qu.:24.75       
##  Max.   :92.00   Max.   :16.000   Max.   :48.00   Max.   :82.00       
##  glasgow_final  
##  Min.   :11.00  
##  1st Qu.:15.00  
##  Median :15.00  
##  Mean   :14.47  
##  3rd Qu.:15.00  
##  Max.   :15.00  
## 
## N observaciones válidas: 34

Matriz de dispersión con correlaciones

pairs(vars_modelo,
      main = "Matriz de Dispersión – Variables del Modelo",
      pch = 19, col = "steelblue", cex = 0.7,
      lower.panel = function(x, y, ...) {
        points(x, y, pch = 19, col = "steelblue", cex = 0.7)
        r <- round(cor(x, y, method = "spearman", use = "complete.obs"), 2)
        usr <- par("usr")
        text(mean(usr[1:2]), mean(usr[3:4]),
             paste0("ρ = ", r), cex = 1.1, col = "darkred", font = 2)
      })

La diagonal muestra cada variable; el panel inferior reporta el coeficiente de Spearman (ρ) para cada par. Valores altos de ρ con hacer_total_100 identifican los mejores predictores; correlaciones altas entre predictores advierten posible multicolinealidad.

Prueba de normalidad Shapiro-Wilk

vars_nombres <- c("hacer_total_100", "escolaridad", "edad",
                  "dias_hospitalizacion", "glasgow_final")

lineas <- sapply(vars_nombres, function(v) {
  sw <- shapiro.test(vars_modelo[[v]])
  sprintf("%-25s W = %.4f | p = %.4f | %s",
          v, sw$statistic, sw$p.value,
          ifelse(sw$p.value > 0.05, "Normal ✓", "No normal ✗"))
})
cat(paste(c("=== PRUEBAS DE NORMALIDAD (Shapiro-Wilk) ===", "", lineas),
          collapse = "\n"))
## === PRUEBAS DE NORMALIDAD (Shapiro-Wilk) ===
## 
## hacer_total_100           W = 0.9537 | p = 0.1582 | Normal ✓
## escolaridad               W = 0.9276 | p = 0.0266 | No normal ✗
## edad                      W = 0.9203 | p = 0.0165 | No normal ✗
## dias_hospitalizacion      W = 0.8256 | p = 0.0001 | No normal ✗
## glasgow_final             W = 0.5362 | p = 0.0000 | No normal ✗

Matriz de correlaciones de Spearman

cor_matrix <- cor(vars_modelo, method = "spearman", use = "complete.obs")
cor_acer   <- sort(cor_matrix["hacer_total_100", -1], decreasing = TRUE)

cat(paste0(
  "=== CORRELACIONES DE SPEARMAN ===\n\n",
  "Matriz completa:\n",
  paste(capture.output(print(round(cor_matrix, 3))), collapse = "\n"),
  "\n\nCorrelaciones con hacer_total_100 (ordenadas):\n",
  paste(capture.output(print(round(cor_acer, 3))), collapse = "\n")
))
## === CORRELACIONES DE SPEARMAN ===
## 
## Matriz completa:
##                      hacer_total_100 escolaridad   edad dias_hospitalizacion
## hacer_total_100                1.000       0.645  0.050               -0.444
## escolaridad                    0.645       1.000 -0.174               -0.254
## edad                           0.050      -0.174  1.000                0.068
## dias_hospitalizacion          -0.444      -0.254  0.068                1.000
## glasgow_final                  0.184       0.139 -0.055               -0.375
##                      glasgow_final
## hacer_total_100              0.184
## escolaridad                  0.139
## edad                        -0.055
## dias_hospitalizacion        -0.375
## glasgow_final                1.000
## 
## Correlaciones con hacer_total_100 (ordenadas):
##          escolaridad        glasgow_final                 edad 
##                0.645                0.184                0.050 
## dias_hospitalizacion 
##               -0.444

Con estos resultados se identifica cuál variable tiene mayor relación con el ACE-R, lo que orienta la construcción del modelo simple. Además, revisar las correlaciones entre los propios predictores permite detectar si hay variables muy relacionadas entre sí antes de pasar al modelo múltiple.


3. Modelo de Regresión Lineal Simple

Dado que la escolaridad mostró la correlación más alta con el ACE-R en el análisis anterior, se construye el modelo simple usando esta variable como predictor.

a. Ecuación del modelo

\[\widehat{\text{ACE-R}} = \beta_0 + \beta_1 \cdot \text{Escolaridad} + \varepsilon\]

modelo_simple <- lm(hacer_total_100 ~ escolaridad, data = vars_modelo)
summary(modelo_simple)
## 
## Call:
## lm(formula = hacer_total_100 ~ escolaridad, data = vars_modelo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.274 -12.807   1.799  10.431  22.726 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   35.980      6.225    5.78 2.06e-06 ***
## escolaridad    3.049      0.642    4.75 4.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.29 on 32 degrees of freedom
## Multiple R-squared:  0.4135, Adjusted R-squared:  0.3952 
## F-statistic: 22.56 on 1 and 32 DF,  p-value: 4.108e-05
b0 <- round(coef(modelo_simple)[1], 3)
b1 <- round(coef(modelo_simple)[2], 3)

cat(paste0(
  "Ecuación estimada:\n",
  "  ACE-R = ", b0, ifelse(b1 >= 0, " + ", " - "), abs(b1), " × Escolaridad\n\n",
  "Interpretación:\n",
  "  β₀ = ", b0, " → Puntuación ACE-R esperada con 0 años de escolaridad\n",
  "  β₁ = ", b1, " → Por cada año adicional de escolaridad, el ACE-R ",
  ifelse(b1 > 0, "aumenta", "disminuye"), " en promedio ", abs(b1), " puntos"
))
## Ecuación estimada:
##   ACE-R = 35.98 + 3.049 × Escolaridad
## 
## Interpretación:
##   β₀ = 35.98 → Puntuación ACE-R esperada con 0 años de escolaridad
##   β₁ = 3.049 → Por cada año adicional de escolaridad, el ACE-R aumenta en promedio 3.049 puntos

b. Análisis del modelo

sm       <- summary(modelo_simple)
p_modelo <- pf(sm$fstatistic[1], sm$fstatistic[2],
               sm$fstatistic[3], lower.tail = FALSE)

cat(paste0(
  "=== BONDAD DE AJUSTE ===\n\n",
  "R²          = ", round(sm$r.squared, 4),
  "  →  Variabilidad explicada: ", round(sm$r.squared * 100, 1), "%\n",
  "R² ajustado = ", round(sm$adj.r.squared, 4), "\n",
  "F           = ", round(sm$fstatistic[1], 3),
  "  |  p-valor = ", round(p_modelo, 4),
  "  |  ", ifelse(p_modelo < 0.05,
                  "Modelo globalmente significativo ✓",
                  "Modelo NO significativo ✗"),
  "\n\nCoeficientes:\n",
  paste(capture.output(print(round(coef(summary(modelo_simple)), 4))),
        collapse = "\n")
))
## === BONDAD DE AJUSTE ===
## 
## R²          = 0.4135  →  Variabilidad explicada: 41.3%
## R² ajustado = 0.3952
## F           = 22.559  |  p-valor = 0  |  Modelo globalmente significativo ✓
## 
## Coeficientes:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  35.9797     6.2251  5.7798        0
## escolaridad   3.0491     0.6420  4.7497        0

Dado que se trabaja con datos clínicos donde hay mucha variabilidad entre pacientes, el R² puede no ser muy alto. Lo importante es evaluar si el coeficiente de escolaridad es estadísticamente significativo y verificar que los supuestos del modelo se cumplan.

c. Evaluación — Supuestos de Gauss-Markov

i. Linealidad
plot(vars_modelo$escolaridad, residuals(modelo_simple),
     xlab = "Escolaridad (años)",
     ylab = "Residuos",
     main = "Residuos vs Predictor — Linealidad",
     pch = 19, col = "steelblue", cex = 1.2)
abline(h = 0, col = "red", lwd = 2, lty = 2)
lines(lowess(vars_modelo$escolaridad, residuals(modelo_simple)),
      col = "orange", lwd = 2)

Si los residuos se distribuyen de forma aleatoria alrededor de cero y la línea naranja es aproximadamente plana, podemos asumir que la relación entre las variables es lineal.

ii. Independencia de errores (Durbin-Watson)
dw <- dwtest(modelo_simple)

cat(paste0(
  "=== TEST DE DURBIN-WATSON ===\n\n",
  "Estadístico DW = ", round(dw$statistic, 4), "\n",
  "p-valor        = ", round(dw$p.value, 4), "\n",
  "Interpretación : ", ifelse(dw$p.value > 0.05,
    "No se detecta autocorrelación en los residuos ✓",
    "Posible autocorrelación en los residuos ✗"), "\n\n",
  "Nota: DW ≈ 2 indica independencia; DW < 1.5 sugiere autocorrelación positiva."
))
## === TEST DE DURBIN-WATSON ===
## 
## Estadístico DW = 1.5093
## p-valor        = 0.0739
## Interpretación : No se detecta autocorrelación en los residuos ✓
## 
## Nota: DW ≈ 2 indica independencia; DW < 1.5 sugiere autocorrelación positiva.
iii. Homocedasticidad (Breusch-Pagan)
bp <- bptest(modelo_simple)

cat(paste0(
  "=== TEST DE BREUSCH-PAGAN ===\n\n",
  "BP = ", round(bp$statistic, 4), "  |  p-valor = ", round(bp$p.value, 4), "\n",
  "Interpretación: ", ifelse(bp$p.value > 0.05,
    "Varianza constante de los errores (homocedasticidad) ✓",
    "Heteroscedasticidad detectada ✗")
))
## === TEST DE BREUSCH-PAGAN ===
## 
## BP = 4.6764  |  p-valor = 0.0306
## Interpretación: Heteroscedasticidad detectada ✗
plot(modelo_simple, which = 3,
     main = "Scale-Location (Homocedasticidad)",
     col = "steelblue", pch = 19)

En este gráfico se espera que la línea roja sea aproximadamente plana, lo que indicaría que la varianza de los errores es constante a lo largo de los valores ajustados. Si la línea sube, el modelo no cumpliría el supuesto de homocedasticidad.

iv. Media cero del error
media_res <- mean(residuals(modelo_simple))

cat(paste0(
  "=== MEDIA DE LOS RESIDUOS ===\n\n",
  "Media de residuos: ", formatC(media_res, format = "e", digits = 4), "\n",
  "Interpretación: Por construcción algebraica del MCO con intercepto, ",
  "la media de residuos es exactamente 0 ✓"
))
## === MEDIA DE LOS RESIDUOS ===
## 
## Media de residuos: 3.3919e-16
## Interpretación: Por construcción algebraica del MCO con intercepto, la media de residuos es exactamente 0 ✓
v. Normalidad de residuos
par(mfrow = c(1, 2))

qqnorm(residuals(modelo_simple),
       main = "QQ-Plot de Residuos",
       pch = 19, col = "steelblue", cex = 1.2)
qqline(residuals(modelo_simple), col = "red", lwd = 2)

hist(residuals(modelo_simple),
     main = "Distribución de Residuos",
     xlab = "Residuos", col = "steelblue", border = "white",
     probability = TRUE, breaks = 8)
curve(dnorm(x, mean = 0, sd = sd(residuals(modelo_simple))),
      add = TRUE, col = "red", lwd = 2)

par(mfrow = c(1, 1))

sw_res <- shapiro.test(residuals(modelo_simple))

cat(paste0(
  "=== NORMALIDAD DE RESIDUOS (Shapiro-Wilk) ===\n\n",
  "W = ", round(sw_res$statistic, 4),
  "  |  p-valor = ", round(sw_res$p.value, 4), "\n",
  "Interpretación: ", ifelse(sw_res$p.value > 0.05,
    "Los residuos siguen distribución normal ✓",
    "Los residuos NO siguen distribución normal ✗")
))
## === NORMALIDAD DE RESIDUOS (Shapiro-Wilk) ===
## 
## W = 0.9581  |  p-valor = 0.2134
## Interpretación: Los residuos siguen distribución normal ✓

El QQ-plot permite ver visualmente si los residuos se comportan de forma normal. Si los puntos siguen la línea roja y el histograma tiene forma de campana, se puede asumir que los errores del modelo siguen una distribución normal.


4. Modelo de Regresión Lineal Múltiple

Además de la escolaridad, se incluyen la edad, los días de hospitalización y el Glasgow final, ya que desde el punto de vista clínico estas variables también podrían estar relacionadas con el desempeño cognitivo de los pacientes con TCE.

a. Ecuación del modelo

\[\widehat{\text{ACE-R}} = \beta_0 + \beta_1 \cdot \text{Escol} + \beta_2 \cdot \text{Edad} + \beta_3 \cdot \text{Días\_Hosp} + \beta_4 \cdot \text{Glasgow\_Final} + \varepsilon\]

modelo_multiple <- lm(hacer_total_100 ~ escolaridad + edad +
                        dias_hospitalizacion + glasgow_final,
                      data = vars_modelo)
summary(modelo_multiple)
## 
## Call:
## lm(formula = hacer_total_100 ~ escolaridad + edad + dias_hospitalizacion + 
##     glasgow_final, data = vars_modelo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.706 -10.548  -1.521  10.525  27.781 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -11.5167    34.4805  -0.334    0.741    
## escolaridad            2.9319     0.6120   4.791 4.54e-05 ***
## edad                   0.3629     0.2699   1.345    0.189    
## dias_hospitalizacion  -0.2706     0.1673  -1.617    0.117    
## glasgow_final          2.9926     2.2406   1.336    0.192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.15 on 29 degrees of freedom
## Multiple R-squared:  0.5448, Adjusted R-squared:  0.482 
## F-statistic: 8.676 on 4 and 29 DF,  p-value: 9.853e-05

b. Análisis del modelo — revisión de variables predictoras

sm_m   <- summary(modelo_multiple)
p_mult <- pf(sm_m$fstatistic[1], sm_m$fstatistic[2],
             sm_m$fstatistic[3], lower.tail = FALSE)
cf_tab <- round(coef(summary(modelo_multiple)), 4)
betas  <- coef(modelo_multiple)

dir_clinica <- c(
  escolaridad          = "A mayor escolaridad, mejor desempeño cognitivo esperado",
  edad                 = "La edad podría influir en la recuperación cognitiva tras el TCE",
  dias_hospitalizacion = "Una hospitalización más larga puede reflejar mayor gravedad de la lesión",
  glasgow_final        = "Mayor puntaje de Glasgow al egreso se asocia con mejor función neurológica"
)

interp <- paste(sapply(names(betas)[-1], function(nm) {
  sig <- cf_tab[nm, "Pr(>|t|)"] < 0.05
  sprintf("  β(%s) = %.3f  %s\n  → %s",
          nm, betas[nm],
          ifelse(sig, "[significativo]", "[no significativo]"),
          dir_clinica[nm])
}), collapse = "\n\n")

cat(paste0(
  "=== BONDAD DE AJUSTE — MODELO MÚLTIPLE ===\n\n",
  "R²          = ", round(sm_m$r.squared, 4),
  "  →  Variabilidad explicada: ", round(sm_m$r.squared * 100, 1), "%\n",
  "R² ajustado = ", round(sm_m$adj.r.squared, 4), "\n",
  "F           = ", round(sm_m$fstatistic[1], 3),
  "  |  p-valor = ", round(p_mult, 4),
  "  |  ", ifelse(p_mult < 0.05,
                  "Modelo globalmente significativo ✓",
                  "Modelo NO significativo ✗"),
  "\n\n=== COEFICIENTES ===\n",
  paste(capture.output(print(cf_tab)), collapse = "\n"),
  "\n\n=== INTERPRETACIÓN DE COEFICIENTES ===\n\n",
  interp
))
## === BONDAD DE AJUSTE — MODELO MÚLTIPLE ===
## 
## R²          = 0.5448  →  Variabilidad explicada: 54.5%
## R² ajustado = 0.482
## F           = 8.676  |  p-valor = 1e-04  |  Modelo globalmente significativo ✓
## 
## === COEFICIENTES ===
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)          -11.5167    34.4805 -0.3340   0.7408
## escolaridad            2.9319     0.6120  4.7907   0.0000
## edad                   0.3629     0.2699  1.3447   0.1892
## dias_hospitalizacion  -0.2706     0.1673 -1.6175   0.1166
## glasgow_final          2.9926     2.2406  1.3356   0.1921
## 
## === INTERPRETACIÓN DE COEFICIENTES ===
## 
##   β(escolaridad) = 2.932  [significativo]
##   → A mayor escolaridad, mejor desempeño cognitivo esperado
## 
##   β(edad) = 0.363  [no significativo]
##   → La edad podría influir en la recuperación cognitiva tras el TCE
## 
##   β(dias_hospitalizacion) = -0.271  [no significativo]
##   → Una hospitalización más larga puede reflejar mayor gravedad de la lesión
## 
##   β(glasgow_final) = 2.993  [no significativo]
##   → Mayor puntaje de Glasgow al egreso se asocia con mejor función neurológica
vif_vals <- vif(modelo_multiple)

cat(paste0(
  "=== FACTOR DE INFLACIÓN DE VARIANZA (VIF) ===\n",
  "VIF > 5: multicolinealidad moderada  |  VIF > 10: multicolinealidad severa\n\n",
  paste(capture.output(print(round(vif_vals, 3))), collapse = "\n"),
  "\n\nConclusión: ", ifelse(all(vif_vals < 5),
    "Todos los VIF < 5: no hay multicolinealidad problemática ✓",
    "Al menos un VIF ≥ 5: revisar variables correlacionadas ✗")
))
## === FACTOR DE INFLACIÓN DE VARIANZA (VIF) ===
## VIF > 5: multicolinealidad moderada  |  VIF > 10: multicolinealidad severa
## 
##          escolaridad                 edad dias_hospitalizacion 
##                1.061                1.044                1.092 
##        glasgow_final 
##                1.064 
## 
## Conclusión: Todos los VIF < 5: no hay multicolinealidad problemática ✓

5. Selección del Modelo — Criterio AIC

Para decidir cuál modelo es el más adecuado se usa el Criterio de Información de Akaike (AIC), que permite comparar modelos teniendo en cuenta tanto el ajuste a los datos como el número de variables incluidas. El modelo con el AIC más bajo es el que se prefiere.

\[\text{AIC} = -2\ln(\hat{L}) + 2k\]

m1 <- lm(hacer_total_100 ~ escolaridad,                                      data = vars_modelo)
m2 <- lm(hacer_total_100 ~ escolaridad + edad,                               data = vars_modelo)
m3 <- lm(hacer_total_100 ~ escolaridad + dias_hospitalizacion,               data = vars_modelo)
m4 <- lm(hacer_total_100 ~ escolaridad + glasgow_final,                      data = vars_modelo)
m5 <- lm(hacer_total_100 ~ escolaridad + edad + dias_hospitalizacion,        data = vars_modelo)
m6 <- lm(hacer_total_100 ~ escolaridad + edad + glasgow_final,               data = vars_modelo)
m7 <- lm(hacer_total_100 ~ escolaridad + edad + dias_hospitalizacion + glasgow_final, data = vars_modelo)

modelos <- list(
  "escolaridad"                                = m1,
  "escolaridad + edad"                         = m2,
  "escolaridad + dias_hosp"                    = m3,
  "escolaridad + glasgow"                      = m4,
  "escolaridad + edad + dias_hosp"             = m5,
  "escolaridad + edad + glasgow"               = m6,
  "escolaridad + edad + dias_hosp + glasgow"   = m7
)

tabla_aic <- data.frame(
  Modelo   = names(modelos),
  k        = sapply(modelos, function(m) length(coef(m))),
  AIC      = round(sapply(modelos, AIC), 2),
  BIC      = round(sapply(modelos, BIC), 2),
  R2_aj    = round(sapply(modelos, function(m) summary(m)$adj.r.squared), 4)
)
tabla_aic <- tabla_aic[order(tabla_aic$AIC), ]

mejor <- tabla_aic$Modelo[1]

cat(paste0(
  "=== TABLA COMPARATIVA DE MODELOS (ordenada por AIC) ===\n\n",
  paste(capture.output(print(tabla_aic, row.names = FALSE)), collapse = "\n"),
  "\n\n*** MEJOR MODELO (menor AIC): ", mejor, " ***"
))
## === TABLA COMPARATIVA DE MODELOS (ordenada por AIC) ===
## 
##                                    Modelo k    AIC    BIC  R2_aj
##  escolaridad + edad + dias_hosp + glasgow 5 283.28 292.44 0.4820
##            escolaridad + edad + dias_hosp 4 283.31 290.94 0.4685
##                   escolaridad + dias_hosp 3 283.37 289.48 0.4534
##              escolaridad + edad + glasgow 4 284.21 291.85 0.4541
##                     escolaridad + glasgow 3 284.75 290.85 0.4308
##                        escolaridad + edad 3 285.25 291.36 0.4223
##                               escolaridad 2 285.89 290.47 0.3952
## 
## *** MEJOR MODELO (menor AIC): escolaridad + edad + dias_hosp + glasgow ***
modelo_nulo <- lm(hacer_total_100 ~ 1, data = vars_modelo)
modelo_full <- lm(hacer_total_100 ~ escolaridad + edad +
                    dias_hospitalizacion + glasgow_final, data = vars_modelo)

modelo_step <- step(modelo_nulo,
                    scope = list(lower = modelo_nulo, upper = modelo_full),
                    direction = "both",
                    trace = TRUE)
## Start:  AIC=203.55
## hacer_total_100 ~ 1
## 
##                        Df Sum of Sq     RSS    AIC
## + escolaridad           1    5276.7  7484.9 187.41
## + dias_hospitalizacion  1    1683.1 11078.5 200.74
## + glasgow_final         1    1132.3 11629.2 202.39
## <none>                              12761.6 203.55
## + edad                  1     153.3 12608.3 205.14
## 
## Step:  AIC=187.41
## hacer_total_100 ~ escolaridad
## 
##                        Df Sum of Sq     RSS    AIC
## + dias_hospitalizacion  1     932.3  6552.6 184.88
## + glasgow_final         1     661.4  6823.4 186.26
## + edad                  1     559.9  6925.0 186.76
## <none>                               7484.9 187.41
## - escolaridad           1    5276.7 12761.6 203.55
## 
## Step:  AIC=184.88
## hacer_total_100 ~ escolaridad + dias_hospitalizacion
## 
##                        Df Sum of Sq     RSS    AIC
## + edad                  1     385.9  6166.7 184.82
## + glasgow_final         1     381.1  6171.5 184.85
## <none>                               6552.6 184.88
## - dias_hospitalizacion  1     932.3  7484.9 187.41
## - escolaridad           1    4525.8 11078.5 200.74
## 
## Step:  AIC=184.82
## hacer_total_100 ~ escolaridad + dias_hospitalizacion + edad
## 
##                        Df Sum of Sq     RSS    AIC
## + glasgow_final         1     357.4  5809.3 184.79
## <none>                               6166.7 184.82
## - edad                  1     385.9  6552.6 184.88
## - dias_hospitalizacion  1     758.3  6925.0 186.76
## - escolaridad           1    4850.4 11017.1 202.55
## 
## Step:  AIC=184.79
## hacer_total_100 ~ escolaridad + dias_hospitalizacion + edad + 
##     glasgow_final
## 
##                        Df Sum of Sq     RSS    AIC
## <none>                               5809.3 184.79
## - glasgow_final         1     357.4  6166.7 184.82
## - edad                  1     362.2  6171.5 184.85
## - dias_hospitalizacion  1     524.1  6333.4 185.73
## - escolaridad           1    4597.5 10406.9 202.61
cat(paste0(
  "\n--- Resumen del modelo seleccionado por stepwise ---\n\n",
  paste(capture.output(summary(modelo_step)), collapse = "\n")
))
## 
## --- Resumen del modelo seleccionado por stepwise ---
## 
## 
## Call:
## lm(formula = hacer_total_100 ~ escolaridad + dias_hospitalizacion + 
##     edad + glasgow_final, data = vars_modelo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.706 -10.548  -1.521  10.525  27.781 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -11.5167    34.4805  -0.334    0.741    
## escolaridad            2.9319     0.6120   4.791 4.54e-05 ***
## dias_hospitalizacion  -0.2706     0.1673  -1.617    0.117    
## edad                   0.3629     0.2699   1.345    0.189    
## glasgow_final          2.9926     2.2406   1.336    0.192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.15 on 29 degrees of freedom
## Multiple R-squared:  0.5448, Adjusted R-squared:  0.482 
## F-statistic: 8.676 on 4 and 29 DF,  p-value: 9.853e-05
cf_final <- coef(modelo_step)
vars_sel <- names(cf_final)[-1]
sm_final <- summary(modelo_step)

eq <- paste0("ACE-R = ", round(cf_final[1], 2))
for (nm in vars_sel) {
  eq <- paste0(eq, ifelse(cf_final[nm] >= 0, " + ", " - "),
               abs(round(cf_final[nm], 3)), " × ", nm)
}

cat(paste0(
  "=== MODELO FINAL SELECCIONADO ===\n\n",
  "Variables retenidas : ", paste(vars_sel, collapse = ", "), "\n",
  "Ecuación final      : ", eq, "\n",
  "R² ajustado         = ", round(sm_final$adj.r.squared, 4), "\n",
  "AIC                 = ", round(AIC(modelo_step), 2), "\n\n",
  "Nota importante:\n",
  "El AIC evalúa el ajuste global del modelo, no la significancia individual de\n",
  "cada variable. Que una variable sea retenida por el stepwise no implica que su\n",
  "coeficiente sea estadísticamente significativo.\n\n",
  "Al revisar los p-valores, solo la escolaridad obtuvo p < 0.05, siendo la única\n",
  "variable realmente robusta para explicar el ACE-R en esta muestra. Las demás\n",
  "fueron retenidas porque mejoran el ajuste global del modelo, pero no alcanzan\n",
  "significancia estadística por sí solas.\n\n",
  "Con n=34 y cuatro predictores, el modelo también corre riesgo de sobreajuste\n",
  "(overfitting): se ajusta bien a estos datos pero puede no generalizarse bien\n",
  "a nuevas muestras. Con un mayor tamaño muestral los resultados serían más estables."
))
## === MODELO FINAL SELECCIONADO ===
## 
## Variables retenidas : escolaridad, dias_hospitalizacion, edad, glasgow_final
## Ecuación final      : ACE-R = -11.52 + 2.932 × escolaridad - 0.271 × dias_hospitalizacion + 0.363 × edad + 2.993 × glasgow_final
## R² ajustado         = 0.482
## AIC                 = 283.28
## 
## Nota importante:
## El AIC evalúa el ajuste global del modelo, no la significancia individual de
## cada variable. Que una variable sea retenida por el stepwise no implica que su
## coeficiente sea estadísticamente significativo.
## 
## Al revisar los p-valores, solo la escolaridad obtuvo p < 0.05, siendo la única
## variable realmente robusta para explicar el ACE-R en esta muestra. Las demás
## fueron retenidas porque mejoran el ajuste global del modelo, pero no alcanzan
## significancia estadística por sí solas.
## 
## Con n=34 y cuatro predictores, el modelo también corre riesgo de sobreajuste
## (overfitting): se ajusta bien a estos datos pero puede no generalizarse bien
## a nuevas muestras. Con un mayor tamaño muestral los resultados serían más estables.

5. Modelo de Regresión Logística

En la sección anterior modelamos el ACE-R como variable continua. Ahora damos el siguiente paso clínico: traducir ese puntaje a una decisión dicotómica que es la que realmente usa el equipo de neuropsicología en la práctica: ¿el paciente presenta deterioro cognitivo (sí/no)?

El punto de corte clínico ampliamente reportado para el ACE-R es 87 puntos: por debajo de ese valor se considera que existe deterioro cognitivo probable. Al dicotomizar la variable respuesta, el problema deja de ser “predecir cuántos puntos saca el paciente” y pasa a ser “predecir la probabilidad de que un paciente con TCE caiga en zona de deterioro cognitivo”, lo cual exige un modelo de Regresión Logística.

5.1. Planteamiento del problema e hipótesis

Variable respuesta (Y): deterioro_cognitivo — binaria construida como:

  • 1 (Sí): ACE-R < 88 (deterioro cognitivo probable)
  • 0 (No): ACE-R ≥ 88 (rendimiento dentro de lo esperado)

Variables predictoras (las mismas del modelo lineal, para mantener comparabilidad clínica):

Variable Tipo Descripción
escolaridad Continua (años) Reserva cognitiva
edad Continua (años) Madurez/declive neuronal
dias_hospitalizacion Continua (días) Proxy de gravedad clínica
glasgow_final Continua (3–15) Estado neurológico al egreso

Justificación del modelo logístico: la variable respuesta es binaria (deterioro sí/no), por lo que un modelo lineal no es apropiado (no garantiza predicciones entre 0 y 1, no respeta la naturaleza Bernoulli del fenómeno y viola los supuestos de Gauss-Markov por construcción). La regresión logística modela la probabilidad mediante la función logit:

\[\ln\!\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k\]

Hipótesis:

  • H₀: Ninguna de las variables predictoras se asocia con la probabilidad de presentar deterioro cognitivo. β₁ = β₂ = β₃ = β₄ = 0
  • H₁: Al menos una de las variables predictoras se asocia significativamente con la probabilidad de presentar deterioro cognitivo.

5.2. Diccionario de datos y preparación

# Construir variable binaria a partir del ACE-R
datos_log <- vars_modelo
datos_log$deterioro_cognitivo <- ifelse(datos_log$hacer_total_100 < 88, 1, 0)
datos_log$deterioro_cognitivo <- factor(datos_log$deterioro_cognitivo,
                                        levels = c(0, 1),
                                        labels = c("No", "Sí"))

# Verificar distribución de la variable respuesta
tabla_y <- table(datos_log$deterioro_cognitivo)
prop_y  <- prop.table(tabla_y)

cat(paste0(
  "=== DISTRIBUCIÓN DE LA VARIABLE RESPUESTA ===\n\n",
  "Punto de corte clínico ACE-R: 88 puntos\n",
  "Total observaciones válidas : ", nrow(datos_log), "\n\n",
  "Frecuencias:\n",
  paste(capture.output(print(tabla_y)), collapse = "\n"),
  "\n\nProporciones:\n",
  paste(capture.output(print(round(prop_y, 3))), collapse = "\n"),
  "\n\nNota: Si una de las dos categorías es muy minoritaria (<20%), el modelo\n",
  "logístico tendrá dificultad para predecir esa clase y la sensibilidad/\n",
  "especificidad reflejarán ese desbalance."
))
## === DISTRIBUCIÓN DE LA VARIABLE RESPUESTA ===
## 
## Punto de corte clínico ACE-R: 88 puntos
## Total observaciones válidas : 34
## 
## Frecuencias:
## 
## No Sí 
##  3 31 
## 
## Proporciones:
## 
##    No    Sí 
## 0.088 0.912 
## 
## Nota: Si una de las dos categorías es muy minoritaria (<20%), el modelo
## logístico tendrá dificultad para predecir esa clase y la sensibilidad/
## especificidad reflejarán ese desbalance.

5.3. Análisis exploratorio bivariado (Pre-selección)

Antes de ajustar el modelo, evaluamos si cada predictora continua difiere entre los pacientes con y sin deterioro cognitivo. Como el tamaño muestral es pequeño (n=34) y en la sección anterior ya verificamos no normalidad en varias variables, complementamos la prueba T paramétrica con la prueba no paramétrica de Mann-Whitney (Wilcoxon).

preds <- c("escolaridad", "edad", "dias_hospitalizacion", "glasgow_final")

resultados_biv <- do.call(rbind, lapply(preds, function(v) {
  x <- datos_log[[v]]
  g <- datos_log$deterioro_cognitivo

  media_no <- mean(x[g == "No"], na.rm = TRUE)
  media_si <- mean(x[g == "Sí"], na.rm = TRUE)

  tt  <- t.test(x ~ g)
  wt  <- suppressWarnings(wilcox.test(x ~ g))

  data.frame(
    Variable      = v,
    Media_NoDet   = round(media_no, 2),
    Media_SiDet   = round(media_si, 2),
    Diferencia    = round(media_no - media_si, 2),
    p_tStudent    = round(tt$p.value, 4),
    p_MannWhitney = round(wt$p.value, 4),
    Conclusion    = ifelse(min(tt$p.value, wt$p.value) < 0.05,
                           "Diferencia significativa ✓",
                           "Sin diferencia significativa ✗")
  )
}))

cat(paste0(
  "=== PRUEBAS BIVARIADAS: PREDICTORAS CONTINUAS vs DETERIORO COGNITIVO ===\n\n",
  paste(capture.output(print(resultados_biv, row.names = FALSE)), collapse = "\n"),
  "\n\nLectura: variables con diferencia significativa son candidatas naturales\n",
  "para entrar al modelo. La interpretación se sostiene en evidencia bivariada,\n",
  "no en una decisión algorítmica."
))
## === PRUEBAS BIVARIADAS: PREDICTORAS CONTINUAS vs DETERIORO COGNITIVO ===
## 
##              Variable Media_NoDet Media_SiDet Diferencia p_tStudent
##           escolaridad       12.33        8.45       3.88     0.0743
##                  edad       37.33       28.74       8.59     0.1256
##  dias_hospitalizacion       16.00       20.68      -4.68     0.6271
##         glasgow_final       14.00       14.52      -0.52     0.6596
##  p_MannWhitney                     Conclusion
##         0.1885 Sin diferencia significativa ✗
##         0.1069 Sin diferencia significativa ✗
##         0.5837 Sin diferencia significativa ✗
##         0.6239 Sin diferencia significativa ✗
## 
## Lectura: variables con diferencia significativa son candidatas naturales
## para entrar al modelo. La interpretación se sostiene en evidencia bivariada,
## no en una decisión algorítmica.
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (v in preds) {
  boxplot(datos_log[[v]] ~ datos_log$deterioro_cognitivo,
          xlab = "Deterioro cognitivo (ACE-R<87)",
          ylab = v,
          main = paste("Distribución de", v),
          col = c("steelblue", "coral"))
}

par(mfrow = c(1, 1))

Los boxplots permiten ver de un vistazo si las distribuciones de cada predictora se separan entre los grupos con y sin deterioro. Cajas que casi no se solapan sugieren que esa variable aporta información discriminante.

5.4. Modelo de Regresión Logística Simple

Coherentes con el modelo lineal de la sección anterior (donde la escolaridad fue la única variable robustamente significativa), construimos el modelo logístico simple usando la escolaridad como predictora.

a. Ecuación del modelo

\[\ln\!\left(\frac{p_{\text{deterioro}}}{1-p_{\text{deterioro}}}\right) = \beta_0 + \beta_1 \cdot \text{Escolaridad}\]

modelo_log_simple <- glm(deterioro_cognitivo ~ escolaridad,
                         data = datos_log, family = binomial(link = "logit"))
summary(modelo_log_simple)
## 
## Call:
## glm(formula = deterioro_cognitivo ~ escolaridad, family = binomial(link = "logit"), 
##     data = datos_log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   5.4957     2.6412   2.081   0.0375 *
## escolaridad  -0.2980     0.2129  -1.399   0.1617  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20.294  on 33  degrees of freedom
## Residual deviance: 17.512  on 32  degrees of freedom
## AIC: 21.512
## 
## Number of Fisher Scoring iterations: 6

b. Interpretación inferencial (β y Odds Ratio)

b0_l <- coef(modelo_log_simple)[1]
b1_l <- coef(modelo_log_simple)[2]
or_l <- exp(b1_l)
ic_l <- exp(confint.default(modelo_log_simple))
p_b1 <- coef(summary(modelo_log_simple))["escolaridad", "Pr(>|z|)"]

cat(paste0(
  "=== MODELO LOGÍSTICO SIMPLE: deterioro_cognitivo ~ escolaridad ===\n\n",
  "β₀ (intercepto)   = ", round(b0_l, 4), "\n",
  "β₁ (escolaridad)  = ", round(b1_l, 4),
  "  |  p-valor = ", round(p_b1, 4),
  "  |  ", ifelse(p_b1 < 0.05, "Significativo ✓", "No significativo ✗"),
  "\n\n=== ODDS RATIO ===\n",
  "OR (escolaridad)  = ", round(or_l, 4),
  "  →  IC 95% = [", round(ic_l["escolaridad", 1], 3), ", ",
  round(ic_l["escolaridad", 2], 3), "]\n\n",
  "Interpretación clínica:\n",
  "Por cada año adicional de escolaridad, las odds (apuestas) de presentar\n",
  "deterioro cognitivo se multiplican por ", round(or_l, 3), ".\n",
  "Como OR ", ifelse(or_l < 1, "< 1, la escolaridad actúa como factor PROTECTOR:\n",
                              "> 1, la escolaridad actúa como factor de RIESGO:\n"),
  "más años de educación se asocian con ",
  ifelse(or_l < 1, "menor", "mayor"),
  " probabilidad de deterioro cognitivo.\n\n",
  "Importante: el OR opera sobre las odds, no sobre la probabilidad directa.\n",
  "Para traducirlo a probabilidad se usa la curva logística completa."
))
## === MODELO LOGÍSTICO SIMPLE: deterioro_cognitivo ~ escolaridad ===
## 
## β₀ (intercepto)   = 5.4957
## β₁ (escolaridad)  = -0.298  |  p-valor = 0.1617  |  No significativo ✗
## 
## === ODDS RATIO ===
## OR (escolaridad)  = 0.7423  →  IC 95% = [0.489, 1.127]
## 
## Interpretación clínica:
## Por cada año adicional de escolaridad, las odds (apuestas) de presentar
## deterioro cognitivo se multiplican por 0.742.
## Como OR < 1, la escolaridad actúa como factor PROTECTOR:
## más años de educación se asocian con menor probabilidad de deterioro cognitivo.
## 
## Importante: el OR opera sobre las odds, no sobre la probabilidad directa.
## Para traducirlo a probabilidad se usa la curva logística completa.

5.5. Modelo de Regresión Logística Múltiple

Incorporamos las cuatro variables clínicas para evaluar si, en conjunto, mejoran la capacidad de predecir el deterioro.

a. Ecuación del modelo

\[\ln\!\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 \cdot \text{Escol} + \beta_2 \cdot \text{Edad} + \beta_3 \cdot \text{Días\_Hosp} + \beta_4 \cdot \text{Glasgow\_Final}\]

modelo_log_mult <- glm(deterioro_cognitivo ~ escolaridad + edad +
                         dias_hospitalizacion + glasgow_final,
                       data = datos_log, family = binomial(link = "logit"))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(modelo_log_mult)
## 
## Call:
## glm(formula = deterioro_cognitivo ~ escolaridad + edad + dias_hospitalizacion + 
##     glasgow_final, family = binomial(link = "logit"), data = datos_log)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -16.564     22.759  -0.728    0.467
## escolaridad            -7.114     10.188  -0.698    0.485
## edad                   -2.355      3.094  -0.761    0.446
## dias_hospitalizacion    1.610      2.302   0.699    0.484
## glasgow_final          11.825     15.672   0.755    0.451
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20.2936  on 33  degrees of freedom
## Residual deviance:  5.9684  on 29  degrees of freedom
## AIC: 15.968
## 
## Number of Fisher Scoring iterations: 13

b. Coeficientes, Odds Ratios e interpretación

cf_log    <- coef(summary(modelo_log_mult))
or_log    <- exp(coef(modelo_log_mult))
ic_log    <- exp(confint.default(modelo_log_mult))

tabla_or <- data.frame(
  Variable = rownames(cf_log),
  Beta     = round(cf_log[, "Estimate"], 4),
  p_valor  = round(cf_log[, "Pr(>|z|)"], 4),
  OR       = round(or_log, 4),
  IC_inf   = round(ic_log[, 1], 4),
  IC_sup   = round(ic_log[, 2], 4),
  Signif   = ifelse(cf_log[, "Pr(>|z|)"] < 0.05, "✓", "✗")
)

cat(paste0(
  "=== TABLA DE COEFICIENTES Y ODDS RATIOS — MODELO MÚLTIPLE ===\n\n",
  paste(capture.output(print(tabla_or, row.names = FALSE)), collapse = "\n"),
  "\n\n=== LECTURA DE LOS ODDS RATIOS ===\n\n",
  "OR > 1  →  la variable AUMENTA las odds de deterioro cognitivo (riesgo)\n",
  "OR < 1  →  la variable DISMINUYE las odds de deterioro cognitivo (protección)\n",
  "OR ≈ 1  →  la variable no modifica las odds\n\n",
  "Si el IC 95% del OR contiene el 1, no hay evidencia estadística de efecto."
))
## === TABLA DE COEFICIENTES Y ODDS RATIOS — MODELO MÚLTIPLE ===
## 
##              Variable     Beta p_valor          OR IC_inf       IC_sup Signif
##           (Intercept) -16.5643  0.4667      0.0000 0.0000 1.509078e+12      ✗
##           escolaridad  -7.1141  0.4850      0.0008 0.0000 3.820628e+05      ✗
##                  edad  -2.3553  0.4465      0.0949 0.0002 4.078010e+01      ✗
##  dias_hospitalizacion   1.6098  0.4844      5.0018 0.0549 4.559479e+02      ✗
##         glasgow_final  11.8253  0.4505 136670.6197 0.0000 2.990184e+18      ✗
## 
## === LECTURA DE LOS ODDS RATIOS ===
## 
## OR > 1  →  la variable AUMENTA las odds de deterioro cognitivo (riesgo)
## OR < 1  →  la variable DISMINUYE las odds de deterioro cognitivo (protección)
## OR ≈ 1  →  la variable no modifica las odds
## 
## Si el IC 95% del OR contiene el 1, no hay evidencia estadística de efecto.
direcciones <- c(
  escolaridad          = "Reserva cognitiva: más años de educación → menor riesgo de deterioro",
  edad                 = "Edad: mayor edad podría asociarse con menor recuperación neuronal",
  dias_hospitalizacion = "Estancia prolongada: proxy de mayor gravedad clínica → mayor riesgo",
  glasgow_final        = "Glasgow alto al egreso: mejor estado neurológico → menor riesgo"
)

interp_log <- paste(sapply(names(or_log)[-1], function(nm) {
  sig <- cf_log[nm, "Pr(>|z|)"] < 0.05
  sprintf("  %s: OR = %.3f  %s\n  → %s",
          nm, or_log[nm],
          ifelse(sig, "[significativo p<0.05]", "[no significativo]"),
          direcciones[nm])
}), collapse = "\n\n")

cat(paste0("=== INTERPRETACIÓN CLÍNICA POR VARIABLE ===\n\n", interp_log))
## === INTERPRETACIÓN CLÍNICA POR VARIABLE ===
## 
##   escolaridad: OR = 0.001  [no significativo]
##   → Reserva cognitiva: más años de educación → menor riesgo de deterioro
## 
##   edad: OR = 0.095  [no significativo]
##   → Edad: mayor edad podría asociarse con menor recuperación neuronal
## 
##   dias_hospitalizacion: OR = 5.002  [no significativo]
##   → Estancia prolongada: proxy de mayor gravedad clínica → mayor riesgo
## 
##   glasgow_final: OR = 136670.620  [no significativo]
##   → Glasgow alto al egreso: mejor estado neurológico → menor riesgo

5.6. Selección del modelo — AIC y stepwise

aic_simple <- AIC(modelo_log_simple)
aic_mult   <- AIC(modelo_log_mult)

tabla_aic_log <- data.frame(
  Modelo = c("Simple (solo escolaridad)", "Múltiple (4 predictores)"),
  k      = c(length(coef(modelo_log_simple)), length(coef(modelo_log_mult))),
  AIC    = round(c(aic_simple, aic_mult), 2),
  Deviance = round(c(modelo_log_simple$deviance, modelo_log_mult$deviance), 2)
)

cat(paste0(
  "=== COMPARACIÓN POR AIC ===\n\n",
  paste(capture.output(print(tabla_aic_log, row.names = FALSE)), collapse = "\n"),
  "\n\nMenor AIC = mejor balance entre ajuste y parsimonia."
))
## === COMPARACIÓN POR AIC ===
## 
##                     Modelo k   AIC Deviance
##  Simple (solo escolaridad) 2 21.51    17.51
##   Múltiple (4 predictores) 5 15.97     5.97
## 
## Menor AIC = mejor balance entre ajuste y parsimonia.
modelo_log_step <- step(modelo_log_mult, direction = "both", trace = TRUE)
## Start:  AIC=15.97
## deterioro_cognitivo ~ escolaridad + edad + dias_hospitalizacion + 
##     glasgow_final
## 
##                        Df Deviance    AIC
## <none>                      5.9684 15.968
## - dias_hospitalizacion  1  10.4108 18.411
## - glasgow_final         1  11.3905 19.390
## - escolaridad           1  15.8058 23.806
## - edad                  1  16.0627 24.063
cat(paste0(
  "\n--- Resumen del modelo logístico seleccionado por stepwise ---\n\n",
  paste(capture.output(summary(modelo_log_step)), collapse = "\n")
))
## 
## --- Resumen del modelo logístico seleccionado por stepwise ---
## 
## 
## Call:
## glm(formula = deterioro_cognitivo ~ escolaridad + edad + dias_hospitalizacion + 
##     glasgow_final, family = binomial(link = "logit"), data = datos_log)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -16.564     22.759  -0.728    0.467
## escolaridad            -7.114     10.188  -0.698    0.485
## edad                   -2.355      3.094  -0.761    0.446
## dias_hospitalizacion    1.610      2.302   0.699    0.484
## glasgow_final          11.825     15.672   0.755    0.451
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20.2936  on 33  degrees of freedom
## Residual deviance:  5.9684  on 29  degrees of freedom
## AIC: 15.968
## 
## Number of Fisher Scoring iterations: 13
cf_step  <- coef(summary(modelo_log_step))
vars_step <- rownames(cf_step)[-1]
sig_step  <- vars_step[cf_step[-1, "Pr(>|z|)"] < 0.05]
nosig_step <- setdiff(vars_step, sig_step)

cat(paste0(
  "=== LECTURA CRÍTICA DE LA SELECCIÓN POR STEPWISE ===\n\n",
  "Variables retenidas por el algoritmo : ", paste(vars_step, collapse = ", "), "\n",
  "Variables con p < 0.05 (robustas)    : ",
  ifelse(length(sig_step) == 0, "ninguna", paste(sig_step, collapse = ", ")), "\n",
  "Variables retenidas pero NO signif.  : ",
  ifelse(length(nosig_step) == 0, "ninguna", paste(nosig_step, collapse = ", ")), "\n\n",
  "Advertencia metodológica (idéntica a la del modelo lineal):\n",
  "El AIC mide ajuste global del modelo, no la significancia individual de cada\n",
  "predictor. Que el stepwise retenga una variable NO significa que su coeficiente\n",
  "sea estadísticamente significativo — solo significa que su inclusión mejora el\n",
  "balance ajuste/parsimonia. La verdadera significancia individual se evalúa con\n",
  "los p-valores y los intervalos de confianza de los OR.\n\n",
  "Con n=34 y cuatro predictores, el modelo logístico hereda el mismo riesgo de\n",
  "sobreajuste (overfitting) que el lineal. La regla práctica de '10 eventos por\n",
  "predictor' (regla EPV) sugiere que con esta muestra no deberíamos incluir más\n",
  "de 2–3 predictores con confianza. Los resultados deben leerse como exploratorios\n",
  "y replicarse en una muestra mayor antes de extraer conclusiones clínicas firmes."
))
## === LECTURA CRÍTICA DE LA SELECCIÓN POR STEPWISE ===
## 
## Variables retenidas por el algoritmo : escolaridad, edad, dias_hospitalizacion, glasgow_final
## Variables con p < 0.05 (robustas)    : ninguna
## Variables retenidas pero NO signif.  : escolaridad, edad, dias_hospitalizacion, glasgow_final
## 
## Advertencia metodológica (idéntica a la del modelo lineal):
## El AIC mide ajuste global del modelo, no la significancia individual de cada
## predictor. Que el stepwise retenga una variable NO significa que su coeficiente
## sea estadísticamente significativo — solo significa que su inclusión mejora el
## balance ajuste/parsimonia. La verdadera significancia individual se evalúa con
## los p-valores y los intervalos de confianza de los OR.
## 
## Con n=34 y cuatro predictores, el modelo logístico hereda el mismo riesgo de
## sobreajuste (overfitting) que el lineal. La regla práctica de '10 eventos por
## predictor' (regla EPV) sugiere que con esta muestra no deberíamos incluir más
## de 2–3 predictores con confianza. Los resultados deben leerse como exploratorios
## y replicarse en una muestra mayor antes de extraer conclusiones clínicas firmes.

5.7. Validación: matriz de confusión y métricas de desempeño

Para el modelo logístico, la validación no se hace con los supuestos de Gauss-Markov (esos aplican al modelo lineal), sino evaluando su poder predictivo sobre los datos. Usamos un punto de corte de probabilidad de 0.5 para clasificar.

# Probabilidades predichas por el modelo final (stepwise)
prob_pred <- predict(modelo_log_step, type = "response")

# Clasificación con umbral 0.5
clase_pred <- factor(ifelse(prob_pred >= 0.5, "Sí", "No"),
                     levels = c("No", "Sí"))

# Matriz de confusión (filas = real, columnas = predicho)
mat_conf <- table(Real = datos_log$deterioro_cognitivo,
                  Predicho = clase_pred)

# Métricas
TN <- mat_conf["No", "No"]
FP <- mat_conf["No", "Sí"]
FN <- mat_conf["Sí", "No"]
TP <- mat_conf["Sí", "Sí"]

accuracy    <- (TP + TN) / sum(mat_conf)
sensibilidad <- TP / (TP + FN)   # capacidad de detectar deterioro real
especificidad <- TN / (TN + FP)  # capacidad de descartar deterioro
precision    <- TP / (TP + FP)   # de los que predijo deterioro, cuántos lo eran

cat(paste0(
  "=== MATRIZ DE CONFUSIÓN (umbral = 0.5) ===\n\n",
  paste(capture.output(print(mat_conf)), collapse = "\n"),
  "\n\n=== MÉTRICAS DE DESEMPEÑO ===\n\n",
  sprintf("Exactitud global (Accuracy)   = %.3f  →  %.1f%% de aciertos totales\n",
          accuracy, accuracy * 100),
  sprintf("Sensibilidad (Recall, VP/(VP+FN)) = %.3f  →  detecta el %.1f%% de los casos con deterioro\n",
          sensibilidad, sensibilidad * 100),
  sprintf("Especificidad (VN/(VN+FP))     = %.3f  →  descarta correctamente el %.1f%% de los sanos cognitivos\n",
          especificidad, especificidad * 100),
  sprintf("Precisión (VP/(VP+FP))         = %.3f\n", precision),
  "\n=== LECTURA CLÍNICA ===\n\n",
  "En el contexto del TCE, la SENSIBILIDAD es la métrica más relevante: dejar\n",
  "pasar un caso de deterioro cognitivo (falso negativo) implica privar al paciente\n",
  "de un programa de rehabilitación neuropsicológica que sí necesita. La\n",
  "ESPECIFICIDAD evita derivaciones innecesarias.\n\n",
  "Limitación: las métricas se calcularon sobre los mismos datos con los que se\n",
  "entrenó el modelo (resustitución), por lo que sobreestiman el desempeño real.\n",
  "Una validación cruzada o muestra de prueba independiente daría cifras más bajas\n",
  "y más realistas."
))
## === MATRIZ DE CONFUSIÓN (umbral = 0.5) ===
## 
##     Predicho
## Real No Sí
##   No  2  1
##   Sí  1 30
## 
## === MÉTRICAS DE DESEMPEÑO ===
## 
## Exactitud global (Accuracy)   = 0.941  →  94.1% de aciertos totales
## Sensibilidad (Recall, VP/(VP+FN)) = 0.968  →  detecta el 96.8% de los casos con deterioro
## Especificidad (VN/(VN+FP))     = 0.667  →  descarta correctamente el 66.7% de los sanos cognitivos
## Precisión (VP/(VP+FP))         = 0.968
## 
## === LECTURA CLÍNICA ===
## 
## En el contexto del TCE, la SENSIBILIDAD es la métrica más relevante: dejar
## pasar un caso de deterioro cognitivo (falso negativo) implica privar al paciente
## de un programa de rehabilitación neuropsicológica que sí necesita. La
## ESPECIFICIDAD evita derivaciones innecesarias.
## 
## Limitación: las métricas se calcularon sobre los mismos datos con los que se
## entrenó el modelo (resustitución), por lo que sobreestiman el desempeño real.
## Una validación cruzada o muestra de prueba independiente daría cifras más bajas
## y más realistas.
# Curva logística usando solo escolaridad para visualización
escol_seq <- seq(min(datos_log$escolaridad), max(datos_log$escolaridad),
                 length.out = 200)
pred_curva <- predict(modelo_log_simple,
                      newdata = data.frame(escolaridad = escol_seq),
                      type = "response")

plot(datos_log$escolaridad,
     as.numeric(datos_log$deterioro_cognitivo) - 1,
     pch = 19, col = "steelblue", cex = 1.2,
     xlab = "Años de escolaridad",
     ylab = "Probabilidad de deterioro cognitivo (ACE-R<87)",
     main = "Modelo logístico simple — Escolaridad vs P(deterioro)")
lines(escol_seq, pred_curva, col = "red", lwd = 2)
abline(h = 0.5, lty = 2, col = "gray")
legend("topright",
       legend = c("Datos observados (0/1)", "Curva logística", "Umbral 0.5"),
       col = c("steelblue", "red", "gray"),
       pch = c(19, NA, NA), lty = c(NA, 1, 2), lwd = c(NA, 2, 1),
       cex = 0.85, bg = "white")

La curva en S muestra cómo la probabilidad estimada de deterioro cognitivo cae a medida que aumentan los años de escolaridad. Donde la curva cruza la línea gris (0.5) está el “umbral de escolaridad” por encima del cual el modelo predice ausencia de deterioro.

5.8. Conclusión Estratégica (Toma de Decisiones)

Mensaje para dirección clínica del Hospital Universitario:

El análisis de los 34 pacientes con TCE moderado y severo permite anticipar, con la información disponible al egreso, qué pacientes tienen mayor probabilidad de presentar deterioro cognitivo cuando se les aplica la prueba ACE-R. El hallazgo más sólido y consistente —tanto en el modelo que predice el puntaje como en el que predice el deterioro— es el efecto protector de la escolaridad: los pacientes con más años de educación formal llegan a la evaluación neuropsicológica con menor probabilidad de caer en la zona de deterioro cognitivo. Es la única variable que se sostiene como verdaderamente robusta en ambos análisis.

Las demás variables clínicas evaluadas (edad, días de hospitalización y Glasgow al egreso) muestran tendencias en la dirección clínicamente esperada, pero con esta muestra tan pequeña no podemos afirmar con confianza que sean predictores firmes por sí solos. Aparecen útiles cuando se combinan en el modelo, pero ninguna alcanza el umbral estadístico que exigiríamos para tomar decisiones clínicas basadas exclusivamente en ellas.

Recomendación operativa: identificar tempranamente, ya desde el ingreso hospitalario, a los pacientes con TCE que tienen baja escolaridad (educación primaria o menos) y priorizarlos para programas de rehabilitación neuropsicológica intensiva, sin esperar a que la evaluación cognitiva confirme el deterioro. Esta población es la que tiene la menor “reserva cognitiva” y la que más se beneficia de una intervención preventiva. Para el resto de variables, el modelo aporta una orientación útil pero no debe reemplazar el juicio del neuropsicólogo tratante.

Salvedad técnica para futuras versiones del estudio: la muestra de 34 pacientes limita la fuerza de las conclusiones y abre espacio a sobreajuste; replicar este análisis con un tamaño muestral mayor (idealmente n ≥ 100) permitiría confirmar qué variables clínicas, además de la escolaridad, merecen incorporarse al protocolo de tamizaje.