Este taller ofrece una formación práctica e integral de los principales métodos de dependencia e interdependencia, utilizando R como herramienta de análisis. A través de bases de datos y scripts completamente funcionales, los participantes aprenderán no solo a ejecutar las técnicas, sino también a seleccionar el método más adecuado según la pregunta de investigación, interpretar correctamente los resultados y presentar hallazgos mediante tablas y gráficos con calidad de publicación científica.

En este taller, usaremos la información de una muestra de la encuesta National Health and Nutrition Examination Survey (NHANES) -2017. Esta información consta de variables sociodemográficas, antropométricas, cardiovasculares, de salud mental (PHQ-9) y estilo de vida.

El Patient Health Questionnaire-9 (PHQ-9) es uno de los instrumentos más utilizados para evaluar síntomas depresivos. Está conformado por nueve preguntas, cada una de las cuales evalúa la frecuencia con la que el participante ha experimentado un determinado síntoma durante las dos últimas semanas.

Todos los ítems utilizan la misma escala de respuesta:

• 0: Para nada
• 1: Varios días
• 2: Más de la mitad de los días
• 3: Casi todos los días

A partir de las variables disponibles, formule entre tres y cinco preguntas de investigación que considere podrían responderse utilizando esta base de datos.

Pregunta 1: _______________________________________________________________

Pregunta 2: _______________________________________________________________

Pregunta 3: _______________________________________________________________

Pregunta 4: _______________________________________________________________

Pregunta 5 _______________________________________________________________

Como pudimos ver, no existe una única respuesta correcta. La pregunta Depende del objetivo del estudio, una misma base de datos puede utilizarse para responder diferentes preguntas de investigación. A continuación se presentan algunos ejemplos que serán desarrollados durante el taller.

• Datos crudos
• Datos estandarizados
• Distancias
• Correlaciones
• Similaridades

library(pacman)
pacman::p_load(
  tidyverse, readxl, labelled, gtsummary, tableone, MVN, 
  FactoMineR, factoextra, psych, caret, pROC, MASS, rpart, 
  rpart.plot, nnet, cluster, sjPlot, GGally, plotly, car, 
  rgl, aplpack, ggcorrplot, ggradar, scales, purrr, patchwork, NeuralNetTools)

curso <- read_excel("base_curso.xlsx")
summary(curso)

##  problema_sueno_consulta salud_general           edad        grupo_edad       
##  Length:9254             Length:9254        Min.   : 0.00   Length:9254       
##  Class :character        Class :character   1st Qu.:11.00   Class :character  
##  Mode  :character        Mode  :character   Median :31.00   Mode  :character  
##                                             Mean   :34.33                     
##                                             3rd Qu.:58.00                     
##                                             Max.   :80.00                     
##                                                                               
##      sexo               raza            educacion         estado_civil      
##  Length:9254        Length:9254        Length:9254        Length:9254       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  quintil_pobreza        talla            peso             imc       
##  Length:9254        Min.   : 78.3   Min.   :  3.20   Min.   :12.30  
##  Class :character   1st Qu.:151.4   1st Qu.: 43.10   1st Qu.:20.40  
##  Mode  :character   Median :161.9   Median : 67.75   Median :25.80  
##                     Mean   :156.6   Mean   : 65.14   Mean   :26.58  
##                     3rd Qu.:171.2   3rd Qu.: 85.60   3rd Qu.:31.30  
##                     Max.   :197.7   Max.   :242.60   Max.   :86.20  
##                     NA's   :1238    NA's   :674      NA's   :1249   
##     cintura            pas1            pad1            hba1c      
##  Min.   : 40.00   Min.   : 72.0   Min.   :  0.00   Min.   : 3.80  
##  1st Qu.: 73.90   1st Qu.:106.0   1st Qu.: 60.00   1st Qu.: 5.20  
##  Median : 91.20   Median :118.0   Median : 70.00   Median : 5.50  
##  Mean   : 89.93   Mean   :121.3   Mean   : 67.84   Mean   : 5.77  
##  3rd Qu.:105.30   3rd Qu.:132.0   3rd Qu.: 76.00   3rd Qu.: 5.90  
##  Max.   :169.50   Max.   :228.0   Max.   :136.00   Max.   :16.20  
##  NA's   :1653     NA's   :2952    NA's   :2952     NA's   :3209   
##    colesterol         hdl         consultas_salud_12m actividad_vigorosa
##  Min.   : 76.0   Min.   : 10.00   Min.   :0.000       Min.   :1.000     
##  1st Qu.:151.0   1st Qu.: 43.00   1st Qu.:1.000       1st Qu.:2.000     
##  Median :176.0   Median : 51.00   Median :2.000       Median :2.000     
##  Mean   :179.9   Mean   : 53.39   Mean   :2.302       Mean   :1.755     
##  3rd Qu.:204.0   3rd Qu.: 61.00   3rd Qu.:3.000       3rd Qu.:2.000     
##  Max.   :446.0   Max.   :189.00   Max.   :8.000       Max.   :2.000     
##  NA's   :2516    NA's   :2516     NA's   :25          NA's   :3398      
##  alcohol_12m          fumo100              phq1               phq2          
##  Length:9254        Length:9254        Length:9254        Length:9254       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##      phq3               phq4               phq5               phq6          
##  Length:9254        Length:9254        Length:9254        Length:9254       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##      phq7               phq8               phq9             phq_total     
##  Length:9254        Length:9254        Length:9254        Min.   : 0.000  
##  Class :character   Class :character   Class :character   1st Qu.: 0.000  
##  Mode  :character   Mode  :character   Mode  :character   Median : 2.000  
##                                                           Mean   : 3.241  
##                                                           3rd Qu.: 5.000  
##                                                           Max.   :25.000  
##                                                           NA's   :4186

Un problema que tiene la base completa, es que algunas variables tienen valores perdidos, por lo que para el análisis descriptivo y de regresión logística, se propone trabajar con una base sin valores perdidos.

base <- na.omit(curso)
summary(base)

##  problema_sueno_consulta salud_general           edad        grupo_edad       
##  Length:3510             Length:3510        Min.   :20.00   Length:3510       
##  Class :character        Class :character   1st Qu.:36.00   Class :character  
##  Mode  :character        Mode  :character   Median :52.00   Mode  :character  
##                                             Mean   :50.84                     
##                                             3rd Qu.:64.00                     
##                                             Max.   :80.00                     
##      sexo               raza            educacion         estado_civil      
##  Length:3510        Length:3510        Length:3510        Length:3510       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##  quintil_pobreza        talla            peso             imc       
##  Length:3510        Min.   :139.7   Min.   : 38.20   Min.   :14.80  
##  Class :character   1st Qu.:159.4   1st Qu.: 67.60   1st Qu.:24.80  
##  Mode  :character   Median :166.5   Median : 79.80   Median :28.60  
##                     Mean   :166.9   Mean   : 83.18   Mean   :29.76  
##                     3rd Qu.:174.0   3rd Qu.: 95.40   3rd Qu.:33.50  
##                     Max.   :195.8   Max.   :191.40   Max.   :67.70  
##     cintura           pas1            pad1            hba1c       
##  Min.   : 62.3   Min.   : 86.0   Min.   :  0.00   Min.   : 4.100  
##  1st Qu.: 88.8   1st Qu.:114.0   1st Qu.: 66.00   1st Qu.: 5.300  
##  Median : 99.6   Median :124.0   Median : 72.00   Median : 5.600  
##  Mean   :100.9   Mean   :126.3   Mean   : 72.32   Mean   : 5.835  
##  3rd Qu.:111.4   3rd Qu.:136.0   3rd Qu.: 80.00   3rd Qu.: 6.000  
##  Max.   :169.5   Max.   :224.0   Max.   :124.00   Max.   :16.200  
##    colesterol         hdl         consultas_salud_12m actividad_vigorosa
##  Min.   : 79.0   Min.   : 10.00   Min.   :0.000       Min.   :1.000     
##  1st Qu.:160.0   1st Qu.: 42.00   1st Qu.:1.000       1st Qu.:2.000     
##  Median :185.0   Median : 50.00   Median :2.000       Median :2.000     
##  Mean   :188.6   Mean   : 53.15   Mean   :2.387       Mean   :1.752     
##  3rd Qu.:215.0   3rd Qu.: 62.00   3rd Qu.:3.000       3rd Qu.:2.000     
##  Max.   :446.0   Max.   :178.00   Max.   :8.000       Max.   :2.000     
##  alcohol_12m          fumo100              phq1               phq2          
##  Length:3510        Length:3510        Length:3510        Length:3510       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##      phq3               phq4               phq5               phq6          
##  Length:3510        Length:3510        Length:3510        Length:3510       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##      phq7               phq8               phq9             phq_total     
##  Length:3510        Length:3510        Length:3510        Min.   : 0.000  
##  Class :character   Class :character   Class :character   1st Qu.: 0.000  
##  Mode  :character   Mode  :character   Mode  :character   Median : 2.000  
##                                                           Mean   : 3.178  
##                                                           3rd Qu.: 5.000  
##                                                           Max.   :25.000

var_label(base$problema_sueno_consulta) <- "Problemas de sueño"
var_label(base$salud_general) <- "Estado de salud general"
var_label(base$edad) <- "Edad (años)"
var_label(base$grupo_edad) <- "Grupo de edad"
var_label(base$sexo) <- "Sexo"
var_label(base$raza) <- "Raza/etnia"
var_label(base$educacion) <- "Nivel educativo" 
var_label(base$estado_civil) <- "Estado civil"
var_label(base$quintil_pobreza) <- "Quintil de pobreza"
var_label(base$talla) <- "Talla (cm)"
var_label(base$peso) <- "Peso (kg)"
var_label(base$imc) <- "Índice de masa corporal (kg/m²)"
var_label(base$cintura) <- "Circunferencia de cintura (cm)"
var_label(base$pas1) <- "Presión arterial sistólica (1.ª medición, mmHg)"
var_label(base$pad1) <- "Presión arterial diastólica (1.ª medición, mmHg)"
var_label(base$hba1c) <- "Hemoglobina glucosilada (HbA1c, %)"
var_label(base$colesterol) <- "Colesterol total (mg/dL)"
var_label(base$hdl) <- "Colesterol HDL (mg/dL)"
var_label(base$actividad_vigorosa) <- "Actividad física vigorosa"
var_label(base$alcohol_12m) <- "Consumo de alcohol en los últimos 12 meses"
var_label(base$fumo100) <- "Ha fumado al menos 100 cigarrillos durante la vida"

var_label(base$phq1) <- "Poco interés o placer en hacer las cosas"
var_label(base$phq2) <- "Sentirse decaído, deprimido o sin esperanza"
var_label(base$phq3) <- "Dificultad para dormir o dormir demasiado"
var_label(base$phq4) <- "Sentirse cansado o con poca energía"
var_label(base$phq5) <- "Poco apetito o comer en exceso"
var_label(base$phq6) <- "Sentirse mal consigo mismo o sentirse un fracaso"
var_label(base$phq7) <- "Dificultad para concentrarse"
var_label(base$phq8) <- "Lentitud o inquietud psicomotora"
var_label(base$phq9) <- "Pensamientos de muerte o autolesión"
var_label(base$phq_total) <- "Puntaje total del PHQ-9"

base <- base %>%
  mutate(
    # Problemas de sueño. Esta variable aparece como 1/2
    problema_sueno_consulta = factor(problema_sueno_consulta,levels = c(2, 1), labels = c("No", "Sí")),

    # Sexo. Esta variable ya aparece como Hombre/Mujer
    sexo = factor(sexo, levels = c("Hombre", "Mujer")),

    # Actividad física vigorosa. Esta variable aparece como 1/2
    actividad_vigorosa = factor(actividad_vigorosa, levels = c(1, 2),labels = c("Sí", "No")),

    # Consumo de alcohol. Esta variable ya aparece como "No"/"Sí".
    alcohol_12m = factor(alcohol_12m, levels = c("No", "Sí")),

    # Tabaquismo. Esta variable ya aparece como "No"/"Sí".
    fumo100 = factor(fumo100, levels = c("No", "Sí")),

    # Estado de salud general
      salud_general = factor(salud_general, levels = 1:5,
      labels = c("Excelente", "Muy buena", "Buena", "Regular", "Mala")),
    
    # Raza / Etnia 
    raza = factor(raza,levels = c(1, 2, 3, 4, 6, 7),
                  labels = c("Mex-Amer","Hispano","Blanco","Negro","Asiático","Otra")),
    
      # Nivel educativo 
      educacion = factor(educacion,levels = 1:5,
                         labels = c("EdSecSin9","EdSecSinDip","EdSecDip","EdUnivInc","Univ/Posg")),
    
    # Estado civil
    estado_civil = factor(estado_civil,levels = 1:6,
                          labels = c("Casado/a","Viudo/a","Divorciado/a","Separado/a","Soltero/a","Unión libre"))  )

phq_labels <- c(
  "Nunca",
  "Varios días",
  "Más de la mitad de los días",
  "Casi todos los días"
)

for (i in 1:9) {

  var <- paste0("phq", i)

  base[[var]] <- factor(
    base[[var]],
    levels = 0:3,
    labels = phq_labels,
    ordered = TRUE
  )
}

Para efectos del análisis multivariado, vamos a categorizar algunas de las variables cuantitativas.

base$imccat <- cut(base$imc,breaks = c(-Inf, 18.5, 25, 30, Inf),
                   labels = c("Bajo peso", "Normal", "Sobrepeso", "Obesidad"))

base$hta <- ifelse(base$pas1 >= 140 | base$pad1 >= 90, "Sí", "No")

base$hba1ccat <- cut(base$hba1c, breaks = c(-Inf, 5.7, 6.4, Inf),
                  labels = c("Normal", "Prediabetes", "Diabetes"))

base_modelos <- base %>%
  dplyr::select(
    problema_sueno_consulta,
    edad, grupo_edad,
    sexo, raza, educacion,
    estado_civil, 
    colesterol, hdl,
    quintil_pobreza,
    talla, peso, cintura,
    imc, imccat,
    pas1, pad1, hta,
    hba1c, hba1ccat,
    salud_general,
    actividad_vigorosa,
    fumo100, alcohol_12m,
    phq1, phq2, phq3, phq4, phq5, phq6, phq7, phq8, phq9,
    phq_total
  ) %>%
  na.omit()

De los gráficos de dispersión clásicos a representaciones para más de tres variables de manera simultánea. A continuación presentamos algunos:

• Matriz de Gráficos de Dispersión
• Gráficos 3D Dinámicos / Nubes de Puntos Tridimensionales
• Gráfico de Radar o Estrella:
• Mapas de Calor (Heatmaps) con Dendrogramas
• Gráficos de Burbujas (Bubble Charts)
• Proyecciones Biplot (ACP / ACM)
• Caras de Chernoff
• Superficie de Kriging (Kriging Surface)

Matriz de gráficos bidimensionales que muestra las relaciones pareadas entre todas las variables de un conjunto de datos simultáneamente.

Representaciones en tres ejes (\(X, Y, Z\)) que permiten rotación e interacción para explorar estructuras en tres dimensiones continuas.

Consiste en una técnica de interpolación espacial que permite estimar valores continuos de una variable en un espacio tridimensional a partir de un conjunto de puntos muestreados.

library(gstat)
library(sf)

# 1. Usando las primeras 100 observaciones
datos_sub <- na.omit(base_modelos[1:100, c("edad", "imc", "hba1c")])
datos_sub$id <- 1:nrow(datos_sub) # ID para numerar los puntos

# 2. Convertir a objeto espacial (Edad en X, IMC en Y)
sf_datos <- st_as_sf(datos_sub, coords = c("edad", "imc"))

# 3. Crear grilla regular para interpolación
grid_x <- seq(min(datos_sub$edad) - 5, max(datos_sub$edad) + 5, length.out = 110)
grid_y <- seq(min(datos_sub$imc) - 5, max(datos_sub$imc) + 5, length.out = 60)
grid_df <- expand.grid(edad = grid_x, imc = grid_y)
sf_grid <- st_as_sf(grid_df, coords = c("edad", "imc"))

# 4. Ajustar modelo de Kriging
v_empirico <- variogram(hba1c ~ 1, sf_datos)
v_modelo <- fit.variogram(v_empirico, model = vgm("Sph"))
kriged_res <- krige(hba1c ~ 1, sf_datos, sf_grid, model = v_modelo)

## [using ordinary kriging]

# 5. Convertir la predicción en matriz para contour()
z_matrix <- matrix(kriged_res$var1.pred, nrow = length(grid_x), ncol = length(grid_y))

contour(x = grid_x, y = grid_y, z = z_matrix,  col = "blue",                
        nlevels = 10, xlab = "Edad",  ylab = "IMC",  
        main = "Superficie de Kriging para HbA1c", labcex = 0.8)
points(datos_sub$edad, datos_sub$imc, col = "red", pch = 19, cex = 0.8)
text(datos_sub$edad, datos_sub$imc, labels = datos_sub$id, pos = 4, cex = 0.7, col = "black")

Para detectar patrones globales o grupos: Coordenadas paralelas, Heatmaps con dendrogramas y Biplots.

Para comparar perfiles individuales: Gráficos de radar/estrella y Caras de Chernoff.

Para explorar relaciones directas: Matriz de dispersión y Gráficos de burbujas.

                             TÉCNICAS MULTIVARIADAS
                                       │
           ┌───────────────────────────┴───────────────────────────┐
           ▼                                                       ▼
  TÉCNICAS DE INTERDEPENDENCIA                            TÉCNICAS DE DEPENDENCIA
   (No Supervisadas)                                         (Supervisadas)
           │                                                       │
   ┌───────┴───────┬──────────────┬──────────────┐         ┌───────┴───────┬──────────┬──────────┐
   ▼               ▼              ▼              ▼         ▼               ▼          ▼          ▼
  PCA             AFE          Clustering       ACM     Regresión        LDA        CART        ANN

Para responder preguntas de investigación que buscan explorar relaciones entre variables sin un objetivo de predicción, se utilizan técnicas de interdependencia. Estas técnicas permiten reducir la dimensionalidad de los datos y visualizar patrones complejos.

Por ejemplo, si tenemos un conjunto de variables relacionadas con la calidad de vida, podemos usar técnicas de interdependencia para identificar componentes principales que resuman la información y nos ayuden a entender mejor las relaciones entre las variables o tal vez nos interese evaluar si podemos reducir su dimensionalidad a partir de la creación de nuevas variables.

Es una técnica estadística de reducción de dimensionalidad que transforma un conjunto de variables en un número menor de componentes principales, conservando la mayor parte posible de la variabilidad original. Por ejemplo, permite resumir varias variables relacionadas con la calidad de vida en una escala global.

Generalmente se aplica a variables cuantitativas y se recomienda estandarizarlas o utilizar la matriz de correlaciones cuando tienen escalas diferentes. El ACP identifica componentes a partir de la matriz de covarianzas o correlaciones, ordenándolos según la cantidad de varianza que explican.

Características de los componentes:

1. Son ortogonales entre sí.
2. Se ordenan de mayor a menor varianza explicada.
3. Se busca conservar la mayor variabilidad con el menor número de componentes posible.

Lo que se busca es generar nuevos componentes con la siguiente estructura:

1. Estandarización de Datos
2. Cálculo de la Matriz de Covarianza/Correlación
3. Cálculo de Autovectores y Autovalores
4. Construcción de Componentes Principales
5. Selección del Número de Componentes y Varianza Explicada
6. Carga de los Componentes
7. Visualización de los componentes

Queremos evaluar si podemos reducir estas 5 variables cuantitativas (edad, imc, pas1, hba1c y phq) a un número menor de componentes principales que resuman la información y nos ayuden a entender mejor las relaciones entre ellas.

Preparación de datos

Seleccionamos las variables cuantitativas de interés

basecont <- base_modelos %>%
  dplyr::select(edad, imc, pas1, hba1c, phq_total)

Prueba de Mardia

Calcula los coeficientes de asimetría y curtosis multivariados de Mardia, con su significación estadística. Para normalidad multivariante, ambos valores p (asimetría y curtosis) deben ser superiores a 0.05.

pruebamult <- mvn(basecont, subset = NULL, mvn_test = "mardia")
pruebamult$univariate_normality

##               Test  Variable Statistic p.value    Normality
## 1 Anderson-Darling      edad    34.281  <0.001 ✗ Not normal
## 2 Anderson-Darling       imc    39.789  <0.001 ✗ Not normal
## 3 Anderson-Darling      pas1    38.562  <0.001 ✗ Not normal
## 4 Anderson-Darling     hba1c   284.232  <0.001 ✗ Not normal
## 5 Anderson-Darling phq_total   250.884  <0.001 ✗ Not normal

pruebamult$multivariate_normality

##              Test Statistic p.value     Method          MVN
## 1 Mardia Skewness 10871.141  <0.001 asymptotic ✗ Not normal
## 2 Mardia Kurtosis    88.038  <0.001 asymptotic ✗ Not normal

El ACP es sensible a la escala de las variables (edad, IMC, PAS, HbA1c y el puntaje PHQ no están en las mismas unidades). Se estandarizan explícitamente para dejar claro el paso, aunque PCA(..., scale=TRUE) lo haga internamente más adelante.

basecont_z <- scale(basecont)
head(basecont_z)

##            edad        imc       pas1      hba1c  phq_total
## [1,]  0.2949302 -1.1891559 -0.9561610 -0.1265618 -0.2856061
## [2,]  0.9240602 -0.8798254 -1.1657093 -0.2203409  1.1686650
## [3,]  1.1528347 -1.0204302 -0.7466127  0.3423340  0.1991509
## [4,]  0.5808983  0.1325290 -0.3275162 -0.2203409 -0.2856061
## [5,] -1.6496535 -0.7392206 -0.5370645 -0.6892367  0.4415294
## [6,]  0.5237047  0.8636739  0.3011286 -0.5954575 -0.7703632

eig_decomp <- eigen(mat_cor)
eig_decomp$values   # autovalores

## [1] 1.6894952 1.1211553 0.9132614 0.7900051 0.4860831

eig_decomp$vectors  # autovectores

##            [,1]       [,2]       [,3]        [,4]         [,5]
## [1,] 0.60551967  0.2856812  0.2206271  0.01368339 -0.709132322
## [2,] 0.22128368 -0.6708816 -0.5508129 -0.36074614 -0.259651262
## [3,] 0.58312088  0.1330221  0.1422420 -0.52834771  0.585569016
## [4,] 0.49393093 -0.1811377 -0.2217180  0.76633978  0.294594389
## [5,] 0.01953805 -0.6463733  0.7606163  0.05698268 -0.005970284

acp_FactoMineR <- PCA(basecont, scale.unit = TRUE, ncp = 5, graph = FALSE)
acp_FactoMineR

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 3510 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

La distancia desde el origen a cada proyección se denomina score o puntuación de cada observación.

Criterios habituales:

• Kaiser: retener componentes con autovalor > 1.
• Codo (scree plot): buscar el punto donde la pendiente se aplana.
• Varianza acumulada: retener componentes hasta explicar ~70-80%.

acp_FactoMineR$eig

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1  1.6894952              33.789903                          33.78990
## comp 2  1.1211553              22.423106                          56.21301
## comp 3  0.9132614              18.265228                          74.47824
## comp 4  0.7900051              15.800102                          90.27834
## comp 5  0.4860831               9.721662                         100.00000

fviz_eig(acp_FactoMineR, choice = "eigenvalue", addlabels = TRUE,
         linecolor = "#FC4E07", geom = "line",
         xlab = "Componente Principal", ylab = "Valor Propio")

Las cargas (loadings) indican qué tanto correlaciona cada variable original con cada componente.

acp_FactoMineR <- PCA(basecont, scale.unit = TRUE, ncp = 2, graph = FALSE)
summary(acp_FactoMineR)

## 
## Call:
## PCA(X = basecont, scale.unit = TRUE, ncp = 2, graph = FALSE) 
## 
## 
## Eigenvalues
##                       Dim.1  Dim.2
## Variance              1.689  1.121
## % of var.            33.790 22.423
## Cumulative % of var. 33.790 56.213
## 
## Individuals (the 10 first)
##               Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## 1         |  1.585 | -0.710  0.009  0.201 | -0.963  0.024  0.369 |
## 2         |  2.098 | -0.401  0.003  0.037 |  0.016  0.000  0.000 |
## 3         |  1.757 |  0.210  0.001  0.014 | -0.724  0.013  0.170 |
## 4         |  0.770 |  0.076  0.000  0.010 | -0.258  0.002  0.112 |
## 5         |  2.056 | -1.808  0.055  0.773 |  0.207  0.001  0.010 |
## 6         |  1.435 |  0.375  0.002  0.068 | -0.216  0.001  0.023 |
## 7         |  1.243 | -0.208  0.001  0.028 | -0.987  0.025  0.631 |
## 8         |  1.580 | -0.105  0.000  0.004 | -1.482  0.056  0.880 |
## 9         |  1.415 |  0.990  0.017  0.489 | -0.737  0.014  0.271 |
## 10        |  1.083 | -0.205  0.001  0.036 | -0.251  0.002  0.054 |
## 
## Variables
##              Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## edad      |  0.787 36.665  0.619 | -0.302  8.161  0.092 |
## imc       |  0.288  4.897  0.083 |  0.710 45.008  0.505 |
## pas1      |  0.758 34.003  0.574 | -0.141  1.769  0.020 |
## hba1c     |  0.642 24.397  0.412 |  0.192  3.281  0.037 |
## phq_total |  0.025  0.038  0.001 |  0.684 41.780  0.468 |

Gráfica de variables

fviz_pca_var(acp_FactoMineR, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE)

En la gráfica de variables, los vectores representan las variables originales y su dirección indica cómo contribuyen a los componentes principales. La longitud del vector refleja la importancia de la variable en el componente.

Otro gráfico útil es el de ver cómo se proyectan los individuos en el espacio de los componentes principales, lo que nos permite identificar patrones, agrupamientos o posibles outliers.

fviz_pca_ind(acp_FactoMineR, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)

En este caso, como tenemos tantas observaciones, es recomendable usar un muestreo o seleccionar un subconjunto de individuos para evitar saturar la gráfica.

sub_individuos <- sample(seq_len(nrow(basecont)),size = 100,replace = FALSE)
nombres_submuestra <- rownames(acp_FactoMineR$ind$coord)[sub_individuos]
fviz_pca_ind(acp_FactoMineR,col.ind = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE,  select.ind = list(name = nombres_submuestra))

Finalmente, también en algunos casos puede ser de interés mostrar como se proyectan los individuos en el espacio de los componentes principales, junto con las variables. Esto se conoce como un biplot.

fviz_pca_biplot(acp_FactoMineR,repel = TRUE,col.var = "#2E9FDF",col.ind = "#00AFBB",
                select.ind = list(name = nombres_submuestra))

El Análisis de Correspondencias Múltiples (ACM) es una técnica de estadística multivariante de tipo exploratorio y descriptivo que permite estudiar simultáneamente varias variables categóricas. Su objetivo es representar los datos en un espacio de pocas dimensiones para identificar asociaciones entre las modalidades de las variables y reconocer individuos con perfiles similares.

El ACM puede considerarse una extensión del Análisis de Correspondencias Simple (ACS). Mientras que el ACS estudia la relación entre dos variables categóricas, el ACM permite analizar simultáneamente más de dos variables categóricas observadas sobre los mismos individuos.

• Estudiar asociaciones entre varias variables categóricas.
• Identificar patrones entre las modalidades de las variables.
• Reconocer individuos con perfiles de respuesta similares.
• Reducir la dimensionalidad de los datos.
• Representar las principales asociaciones mediante mapas factoriales.
• Obtener dimensiones que puedan utilizarse posteriormente en otros análisis, como la clasificación.

Matriz disyuntiva completa. Con \(n\) individuos y \(Q\) variables categóricas (la variable \(q\) con \(J_q\) modalidades, \(J=\sum_{q=1}^{Q}J_q\)), se construye una matriz \(\mathbf{X}\) (\(n \times J\)) donde \(x_{ij}=1\) si el individuo \(i\) presenta la modalidad \(j\), y \(0\) en caso contrario. Cada fila suma exactamente \(Q\) unos.

Tabla de Burt. \(\mathbf{B} = \mathbf{X}^{\top}\mathbf{X}\), matriz cuadrada y simétrica (\(J \times J\)). Los bloques fuera de la diagonal recogen asociaciones entre variables; los de la diagonal, las frecuencias de las modalidades.

• Distancia \(\chi^2\): compara perfiles ponderando por frecuencia; las modalidades poco frecuentes pesan más y pueden aparecer alejadas del origen sin ser necesariamente las más importantes.
• Inercia total: \(\dfrac{J-Q}{Q}\); mide la dispersión de los perfiles, de forma análoga a la varianza en ACP.
• Valores propios: \(\lambda_1 \geq \lambda_2 \geq \cdots\), obtenidos por descomposición en valores singulares. Cada uno indica la inercia de su dimensión; el primer eje concentra la mayor parte.
• Correcciones de Benzécri y Greenacre: ajustan los porcentajes de inercia explicada, que en bruto suelen parecer bajos. Al reportarlos, indicar si son originales o corregidos.

• Coordenadas factoriales: posición de individuos/modalidades en los ejes; no se interpretan solas, sino junto con contribución y \(\cos^2\).
• Contribuciones: cuánto influye cada punto en la orientación de un eje.
• Coseno cuadrado (\(\cos^2\)): calidad de representación (cerca de 1 = buena; cerca de 0 = mala).

• Individuos cercanos → perfiles similares.
• Modalidades cercanas → asociadas entre sí.
• Puntos alejados del origen → perfiles diferenciados del promedio (o modalidades poco frecuentes).

1. Revisar la inercia de los ejes.
2. Identificar modalidades con mayor contribución.
3. Comprobar su \(\cos^2\).
4. Observar extremos opuestos de cada eje.
5. Examinar modalidades e individuos próximos entre sí.
6. Complementar con variables suplementarias (no participan en la construcción de los ejes, solo se proyectan después).

• Agrupar categorías de muy baja frecuencia si es necesario.
• Interpretar los ejes según contribución y \(\cos^2\), no solo distancia al origen.
• El ACM es exploratorio/descriptivo, no un contraste de hipótesis.
• Sus coordenadas pueden usarse luego en clasificación (p. ej., HCPC).

Para el ejercicio no vamos a incluir los 9 ítems del PHQ-9 junto con todas las demás variables, ya que tendríamos demasiadas categorías.

base$hba1ccat <- cut(base$hba1c, breaks = c(-Inf, 5.6, 6.4, Inf),
  labels = c("HbA1c_N", "Prediab", "Diab"))

base$imccat <- cut(base$imc, breaks = c(-Inf, 18.5, 24.9, 29.9, Inf), 
                labels = c("B_Peso", "IMC_N", "SOBREP", "OBES"))

base$educacion2 <- factor(base$educacion, labels = c("EdSecSin9", "EdSecSinDip", "EdSecDip", "EdUnivInc", "Univ/Posg")) 

base$raza2 <-  factor(base$raza, labels = c("Mex_Am", "Hisp", "Blanco", "Negro", "Asiático","Otra_raza"))

base$sexo <- factor(base$sexo, labels = c("Masc","Fem"))

base$fuma <- factor(base$fumo100, labels = c("Fum100","NoFum100"))

base$hta <- factor(base$hta, labels = c("HtaNo","HtaSi"))

base$salud_general2 <- factor(base$salud_general, labels=c("SG_Exc","SG_MB","SG_B","SG_Reg","SG_Mala")) 

base$quintil_pobreza2 <- factor(base$quintil_pobreza, 
                                labels=c("QP_M_Alto","QP_Alto","QP_Med","QP_Bajo","QP_M_Bajo"))

base$actividad_vigorosa2 <- factor(base$actividad_vigorosa, labels=c("Act_Vig_Sí","Act_Vig_No"))

base$problema_sueno_consulta2 <- factor(base$problema_sueno_consulta, 
                                        labels=c("Prob_Sueño_Sí", "Prob_Sueño_No"))

library(FactoMineR)
library(factoextra)
library(dplyr)

datos_acm <- base %>%
  dplyr::select(hba1ccat, hta, imccat, educacion2, sexo, raza2, fuma,
         salud_general2, quintil_pobreza2, actividad_vigorosa2,
         problema_sueno_consulta2
  ) %>%
  na.omit()

acm <- MCA(datos_acm, quali.sup = "problema_sueno_consulta2", graph = FALSE)

summary(acm)

## 
## Call:
## MCA(X = datos_acm, quali.sup = "problema_sueno_consulta2", graph = FALSE) 
## 
## 
## Eigenvalues
##                       Dim.1  Dim.2  Dim.3  Dim.4  Dim.5
## Variance              0.229  0.149  0.137  0.125  0.119
## % of var.             8.795  5.740  5.279  4.808  4.574
## Cumulative % of var.  8.795 14.535 19.814 24.622 29.196
## 
## Individuals (the 10 first)
##                  Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## 1             | -1.223  0.186  0.484 |  0.500  0.048  0.081 | -0.048  0.000
## 2             |  0.027  0.000  0.000 | -0.537  0.055  0.170 |  0.357  0.026
## 3             |  0.205  0.005  0.010 | -0.276  0.015  0.019 |  0.243  0.012
## 4             | -0.970  0.117  0.334 |  0.252  0.012  0.023 |  0.077  0.001
## 5             |  0.016  0.000  0.000 | -0.501  0.048  0.136 |  0.485  0.049
## 6             |  0.496  0.031  0.073 |  1.031  0.203  0.317 |  0.349  0.025
## 7             | -0.276  0.009  0.033 | -0.271  0.014  0.032 |  0.250  0.013
## 8             | -0.957  0.114  0.365 |  0.034  0.000  0.000 |  0.499  0.052
## 9             |  0.841  0.088  0.195 |  0.632  0.076  0.110 | -0.215  0.010
## 10            | -0.519  0.034  0.158 | -0.046  0.000  0.001 | -0.368  0.028
##                 cos2  
## 1              0.001 |
## 2              0.075 |
## 3              0.015 |
## 4              0.002 |
## 5              0.128 |
## 6              0.036 |
## 7              0.027 |
## 8              0.099 |
## 9              0.013 |
## 10             0.079 |
## 
## Categories (the 10 first)
##                   Dim.1     ctr    cos2  v.test     Dim.2     ctr    cos2
## HbA1c_N       |  -0.279   1.901   0.098 -18.554 |  -0.138   0.715   0.024
## Prediab       |   0.182   0.438   0.014   7.106 |   0.023   0.011   0.000
## Diab          |   0.722   3.171   0.084  17.192 |   0.505   2.371   0.041
## HtaNo         |  -0.137   0.622   0.059 -14.384 |  -0.031   0.049   0.003
## HtaSi         |   0.431   1.956   0.059  14.384 |   0.097   0.153   0.003
## B_Peso        |   0.142   0.014   0.000   1.051 |  -1.107   1.262   0.019
## IMC_N         |  -0.520   2.853   0.086 -17.367 |  -0.083   0.110   0.002
## SOBREP        |  -0.149   0.313   0.011  -6.087 |   0.363   2.834   0.062
## OBES          |   0.405   3.033   0.120  20.529 |  -0.188   1.004   0.026
## EdSecSin9     |   1.204   4.062   0.099  18.662 |   2.390  24.537   0.391
##                v.test     Dim.3     ctr    cos2  v.test  
## HbA1c_N        -9.196 |   0.564  12.930   0.401  37.492 |
## Prediab         0.897 |  -0.681  10.273   0.203 -26.662 |
## Diab           12.010 |  -0.772   6.041   0.096 -18.383 |
## HtaNo          -3.248 |   0.236   3.080   0.175  24.795 |
## HtaSi           3.248 |  -0.742   9.685   0.175 -24.795 |
## B_Peso         -8.194 |   2.431   6.624   0.092  18.001 |
## IMC_N          -2.760 |   0.638   7.154   0.129  21.307 |
## SOBREP         14.785 |  -0.196   0.894   0.018  -7.965 |
## OBES           -9.544 |  -0.304   2.846   0.068 -15.406 |
## EdSecSin9      37.054 |   0.425   0.843   0.012   6.585 |
## 
## Categorical variables (eta2)
##                       Dim.1 Dim.2 Dim.3  
## hba1ccat            | 0.126 0.046 0.401 |
## hta                 | 0.059 0.003 0.175 |
## imccat              | 0.142 0.078 0.240 |
## educacion2          | 0.511 0.518 0.111 |
## sexo                | 0.002 0.000 0.012 |
## raza2               | 0.220 0.548 0.116 |
## fuma                | 0.136 0.117 0.006 |
## salud_general2      | 0.463 0.119 0.096 |
## quintil_pobreza2    | 0.377 0.054 0.180 |
## actividad_vigorosa2 | 0.251 0.010 0.034 |
## 
## Supplementary categories
##                  Dim.1   cos2 v.test    Dim.2   cos2 v.test    Dim.3   cos2
## Prob_Sueño_Sí | -0.100  0.024 -9.224 |  0.063  0.010  5.816 |  0.025  0.002
## Prob_Sueño_No |  0.241  0.024  9.224 | -0.152  0.010 -5.816 | -0.060  0.002
##               v.test  
## Prob_Sueño_Sí  2.310 |
## Prob_Sueño_No -2.310 |
## 
## Supplementary categorical variables (eta2)
##                            Dim.1 Dim.2 Dim.3  
## problema_sueno_consulta2 | 0.024 0.010 0.002 |

fviz_screeplot(
  acm,
  addlabels = TRUE,
  main = "Varianza explicada por dimensión (ACM)"
)

fviz_mca_var(acm, choice = "mca.cor",
  repel = TRUE, ggtheme = theme_minimal(),
  title = "Correlación entre variables activas y las dimensiones")

fviz_contrib(acm, choice = "var", axes = 1, top = 15) +
  labs(title = "Contribución de las categorías a la Dimensión 1")

fviz_contrib(acm, choice = "var", axes = 2, top = 15) +
  labs(title = "Contribución de las categorías a la Dimensión 2")

fviz_mca_var(acm, repel = TRUE, ggtheme = theme_minimal(),
  title = "ACM - categorías ", cex = 0.7)

fviz_mca_var(acm, select.var = list(contrib = 20),   
  col.var = "contrib", gradient.cols = c("#2C728E", "#3CBB75", "#FDE725"),
  repel = TRUE, labelsize = 3,
  ggtheme = theme_minimal(),
  title = "ACM - Categorías por contribución")

A la izquierda se agrupan personas con mayor nivel educativo, excelente percepción de salud, alta actividad física y menor pobreza, asociadas principalmente con la etnia asiática. A la derecha se concentra un perfil de mayor vulnerabilidad y morbilidad, caracterizado por menor escolaridad, salud regular, diabetes y origen hispano o mexicano-americano, destacando la secundaria incompleta como categoría de mayor contribución.

El análisis factorial (AF) es una técnica estadística multivariada que permite identificar estructuras latentes que explican las correlaciones entre un conjunto de variables observadas. Su principal objetivo es reducir la dimensionalidad e identificar, si existen, factores subyacentes. Esta técnica es ampliamente utilizada en el desarrollo y validación de instrumentos psicométricos, como escalas de ansiedad, calidad de vida, entre otros. El AF se divide en dos grandes enfoques:

El EFA se utiliza cuando no se cuenta con una hipótesis previa sobre la estructura factorial del instrumento. Su objetivo es explorar los datos para descubrir patrones o dimensiones latentes, para lo cual busca identificar cuántos factores existen y qué variables o ítems se asocian a cada uno.

Es de uso común en las etapas iniciales de diseño de cuestionarios o escalas, pero también es ampliamente usado cuando se busca reducir la dimensionalidad de constructos, identificando ítems con cargas débiles o irrelevantes para el mismo.

Nota: No confundir con el Análisis de Componentes Principales (ACP).

El CFA, por el contrario, se aplica cuando existe una hipótesis teórica o empírica sobre la estructura factorial. Esta técnica permite contrastar si los factores o si las cargas factoriales o ambas, se ajustan a un modelo factorial previamente especificado. Es fundamental en la validación estructural de instrumentos y en la evaluación de modelos teóricos.

Modelo matemático

Antes de ejecutar cualquier paso práctico conviene tener claro qué está estimando el EFA. Sea \(X_i\) una variable observada. El modelo factorial se expresa como:

\[ X_i = \lambda_{i1}F_1 + \lambda_{i2}F_2 + \cdots + \lambda_{im}F_m + \varepsilon_i \]

donde:

• \(X_i\): i-ésima variable observada, - \(F_j\): j-ésimo factor común (latente)
• \(\lambda_{ij}\): carga factorial de la variable \(X_i\) en el factor \(F_j\)
• \(\varepsilon_i\): término de error específico de \(X_i\)
• \(m\): número de factores comunes
• \(i = 1, \dots, p\), siendo \(p\) el número de variables observadas

En notación matricial, para todas las variables a la vez:

\[ \mathbf{X} = \boldsymbol{\Lambda} \mathbf{F} + \boldsymbol{\varepsilon} \]

donde \(\mathbf{X}\) es el vector \(p \times 1\) de variables observadas, \(\boldsymbol{\Lambda}\) es la matriz \(p \times m\) de cargas factoriales, \(\mathbf{F}\) es el vector \(m \times 1\) de factores comunes y \(\boldsymbol{\varepsilon}\) es el vector \(p \times 1\) de errores específicos.

\[ \text{Cov}(\mathbf{F}) = \mathbf{I} \]

Los factores latentes son estandarizados y no correlacionados entre sí (varianza = 1, covarianza = 0), lo que facilita la interpretación y la identificación del modelo.

\[ \operatorname{Cov}(\boldsymbol{\varepsilon}) = \boldsymbol{\Psi} \]

\(\boldsymbol{\Psi}\) es diagonal: los errores no están correlacionados entre sí, y cada ítem tiene su propia varianza específica.

\[ \operatorname{Cov}(\mathbf{F}, \boldsymbol{\varepsilon}) = \mathbf{0} \]

Los factores comunes y los errores no están correlacionados, lo que garantiza que la varianza explicada por los factores no se contamine con los errores.

\[ \mathbf{\Sigma} = \boldsymbol{\Lambda} \boldsymbol{\Lambda}^\top + \boldsymbol{\Psi} \]

La diagonal de \(\Lambda\Lambda^\top\) para la variable \(i\) es la comunalidad \((h_i^2 = \sum_j \lambda_{ij}^2)\) y el elemento correspondiente de \(\Psi\) es la unicidad \((\psi_i)\), de modo que \(\text{Var}(X_i) = h_i^2 + \psi_i\): la varianza de cada ítem se descompone en una parte explicada por los factores comunes y otra atribuible a varianza específica o error.

Con el modelo ya claro, el EFA se ejecuta siguiendo estos pasos, que desarrollaremos uno a uno aplicándolos a los nueve ítems del PHQ-9:

1. Confiabilidad preliminar de la escala (Alfa de Cronbach).
2. Evaluar la factibilidad del análisis (tamaño muestral, correlaciones, determinante, KMO, Bartlett).
3. Elegir el método de extracción (PAF, ML, ULS, GLS, regularizado, bayesiano).
4. Determinar el número de factores a retener (Kaiser, scree plot, análisis paralelo, MAP, EGA).
5. Extraer la solución factorial inicial e interpretar sus índices de ajuste.
6. Rotar los factores para mejorar la interpretabilidad (ortogonal vs. oblicua).
7. Interpretar las cargas factoriales resultantes.
8. Refinar el instrumento (eliminar ítems problemáticos, repetir el proceso).

Paso 0: Confiabilidad preliminar de la escala

Inicialmente revisamos la consistencia interna del instrumento con el alfa de Cronbach, para tener una idea de si los ítems están midiendo algo en común. Un valor ≥ 0.70 suele considerarse aceptable.

library(psych)
library(dplyr)

datos_phq <- base %>%
  dplyr::select(phq1, phq2, phq3, phq4, phq5, phq6, phq7, phq8, phq9) %>%
  mutate(across(everything(), ~ as.numeric(.x) - 1))

cat("N usado en el AFE:", nrow(datos_phq), "\n")

## N usado en el AFE: 3510

psych::alpha(datos_phq)

## 
## Reliability analysis   
## Call: psych::alpha(x = datos_phq)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
##       0.82      0.83    0.83      0.36   5 0.0042 0.35 0.46     0.35
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.82  0.82  0.83
## Duhachek  0.82  0.82  0.83
## 
##  Reliability if an item is dropped:
##      raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
## phq1      0.80      0.81    0.81      0.35 4.4   0.0047 0.0082  0.36
## phq2      0.79      0.80    0.79      0.33 4.0   0.0050 0.0058  0.35
## phq3      0.81      0.82    0.81      0.36 4.5   0.0046 0.0078  0.36
## phq4      0.80      0.81    0.80      0.35 4.3   0.0048 0.0075  0.35
## phq5      0.81      0.82    0.81      0.36 4.5   0.0046 0.0086  0.35
## phq6      0.80      0.81    0.80      0.35 4.2   0.0046 0.0077  0.35
## phq7      0.81      0.81    0.81      0.35 4.4   0.0046 0.0092  0.35
## phq8      0.81      0.82    0.82      0.37 4.6   0.0044 0.0092  0.37
## phq9      0.83      0.83    0.82      0.38 4.9   0.0043 0.0053  0.36
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop  mean   sd
## phq1 3510  0.67  0.66  0.60   0.55 0.374 0.75
## phq2 3510  0.76  0.77  0.76   0.67 0.343 0.70
## phq3 3510  0.69  0.64  0.57   0.54 0.632 0.95
## phq4 3510  0.73  0.68  0.63   0.59 0.749 0.92
## phq5 3510  0.66  0.64  0.57   0.53 0.380 0.77
## phq6 3510  0.67  0.70  0.66   0.57 0.232 0.60
## phq7 3510  0.65  0.66  0.59   0.54 0.260 0.67
## phq8 3510  0.56  0.60  0.52   0.47 0.157 0.52
## phq9 3510  0.44  0.53  0.44   0.38 0.051 0.29
## 
## Non missing response frequency for each item
##         0    1    2    3 miss
## phq1 0.75 0.16 0.05 0.04    0
## phq2 0.76 0.17 0.04 0.03    0
## phq3 0.61 0.23 0.07 0.09    0
## phq4 0.50 0.32 0.09 0.08    0
## phq5 0.75 0.15 0.05 0.04    0
## phq6 0.84 0.11 0.03 0.02    0
## phq7 0.83 0.10 0.03 0.03    0
## phq8 0.90 0.06 0.02 0.02    0
## phq9 0.96 0.03 0.01 0.00    0

Tamaño de muestra recomendado

• Se sugiere entre 5 y 10 sujetos por ítem.
• En general, se recomienda un mínimo de 100 a 200 casos.
• Estudios más exigentes indican que una muestra de 300 es buena, 500 muy buena y 1000 excelente.

Correlaciones entre los ítems: se espera que existan correlaciones significativas y moderadas (idealmente ≥ 0.30).

Determinante de la matriz de correlaciones: debe ser > 0.00001; valores cercanos o más pequeños indican multicolinealidad.

KMO (Kaiser-Meyer-Olkin): mide si el número de observaciones y las correlaciones entre variables son suficientes para un análisis factorial confiable.

Valores muy cercanos a 1 también pueden ser síntoma de redundancia

Prueba de esfericidad de Bartlett: contrasta

\[ H_0: R = I \qquad H_1: R \neq I \]

Se espera un p < 0.05, es decir, evidencia de correlación al menos moderada entre los ítems.

matriz_cor <- cor(datos_phq)
KMO(matriz_cor)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = matriz_cor)
## Overall MSA =  0.89
## MSA for each item = 
## phq1 phq2 phq3 phq4 phq5 phq6 phq7 phq8 phq9 
## 0.90 0.86 0.88 0.87 0.92 0.87 0.92 0.92 0.89

cortest.bartlett(matriz_cor, n = nrow(datos_phq))

## $chisq
## [1] 8940.259
## 
## $p.value
## [1] 0
## 
## $df
## [1] 36

Depende del objetivo del análisis y de las características de los datos. Entre las principales alternativas se encuentran la factorización de ejes principales (PAF), máxima verosimilitud (ML), mínimos cuadrados no ponderados (ULS), mínimos cuadrados generalizados (GLS), métodos regularizados y métodos bayesianos.

La PAF es adecuada para identificar estructuras latentes y no requiere normalidad multivariada estricta, ya que trabaja con la varianza común. ML permite realizar pruebas e índices de ajuste, pero requiere asumir normalidad multivariada y es especialmente útil cuando posteriormente se realizará un CFA. ULS es una alternativa robusta frente a distribuciones asimétricas y puede ser apropiada para escalas tipo Likert, mientras que GLS es similar a ML, aunque algo más tolerante frente a desviaciones leves de la normalidad. Los métodos regularizados, como LASSO, penalizan cargas pequeñas para obtener estructuras más simples, y los métodos bayesianos incorporan información previa y pueden ser útiles en muestras pequeñas o cuando existen problemas de convergencia.

El análisis de factores comunes incluye métodos como PAF, ML, ULS y GLS, y se diferencia del PCA porque busca explicar principalmente la varianza común entre las variables, mientras que el PCA utiliza la varianza total.

Para los ítems del PHQ-9, que son variables ordinales tipo Likert y no cumplen necesariamente una normalidad multivariada estricta, se opta por la factorización de ejes principales (PAF; fm = "pa").

Debe determinarse combinando varios criterios. Kaiser es sencillo, pero puede sobreestimar; el scree plot identifica el “codo”, pero es subjetivo. El análisis paralelo de Horn y la prueba MAP de Velicer ofrecen criterios más robustos para evitar la sobreextracción. EGANet puede complementar el análisis en estructuras complejas. La decisión final debe integrar evidencia estadística, interpretabilidad y teória.

resultado_paralelo <- fa.parallel(datos_phq, fa = "fa", fm = "pa")

## Parallel analysis suggests that the number of factors =  3  and the number of components =  NA

nfactores_sugeridos <- resultado_paralelo$nfact
cat("Número de factores sugerido por el análisis paralelo:", nfactores_sugeridos, "\n")

## Número de factores sugerido por el análisis paralelo: 3

modelo_afe <- fa(matriz_cor, nfactors = 3, fm = "pa", rotate = "none")
salida <- capture.output(print(modelo_afe, cut = 0, digits = 3))
salida

##  [1] "Factor Analysis using method =  pa"                                    
##  [2] "Call: fa(r = matriz_cor, nfactors = 3, rotate = \"none\", fm = \"pa\")"
##  [3] "Standardized loadings (pattern matrix) based upon correlation matrix"  
##  [4] "       PA1    PA2    PA3    h2    u2  com"                             
##  [5] "phq1 0.614 -0.030 -0.194 0.415 0.585 1.20"                             
##  [6] "phq2 0.786  0.195 -0.213 0.702 0.298 1.28"                             
##  [7] "phq3 0.586 -0.317  0.067 0.449 0.551 1.57"                             
##  [8] "phq4 0.652 -0.366 -0.025 0.560 0.440 1.58"                             
##  [9] "phq5 0.571 -0.172  0.045 0.357 0.643 1.19"                             
## [10] "phq6 0.672  0.271 -0.008 0.524 0.476 1.32"                             
## [11] "phq7 0.594  0.047  0.122 0.370 0.630 1.10"                             
## [12] "phq8 0.527  0.104  0.244 0.348 0.652 1.50"                             
## [13] "phq9 0.441  0.281  0.094 0.282 0.718 1.81"                             
## [14] ""                                                                      
## [15] "                        PA1   PA2   PA3"                               
## [16] "SS loadings           3.364 0.468 0.174"                               
## [17] "Proportion Var        0.374 0.052 0.019"                               
## [18] "Cumulative Var        0.374 0.426 0.445"                               
## [19] "Proportion Explained  0.840 0.117 0.043"                               
## [20] "Cumulative Proportion 0.840 0.957 1.000"                               
## [21] ""                                                                      
## [22] "Mean item complexity =  1.4"                                           
## [23] "Test of the hypothesis that 3 factors are sufficient."                 
## [24] ""                                                                      
## [25] "df null model =  36  with the objective function =  2.551"             
## [26] "df of  the model are 12  and the objective function was  0.02 "        
## [27] ""                                                                      
## [28] "The root mean square of the residuals (RMSR) is  0.014 "               
## [29] "The df corrected root mean square of the residuals is  0.024 "         
## [30] ""                                                                      
## [31] "Fit based upon off diagonal values = 0.999"                            
## [32] "Measures of factor score adequacy             "                        
## [33] "                                                    PA1    PA2    PA3" 
## [34] "Correlation of (regression) scores with factors   0.935  0.696  0.504" 
## [35] "Multiple R square of scores with factors          0.874  0.484  0.254" 
## [36] "Minimum correlation of possible factor scores     0.748 -0.032 -0.492"

print(modelo_afe$loadings, cutoff = 0.30)

## 
## Loadings:
##      PA1    PA2    PA3   
## phq1  0.614              
## phq2  0.786              
## phq3  0.586 -0.317       
## phq4  0.652 -0.366       
## phq5  0.571              
## phq6  0.672              
## phq7  0.594              
## phq8  0.527              
## phq9  0.441              
## 
##                  PA1   PA2   PA3
## SS loadings    3.364 0.468 0.174
## Proportion Var 0.374 0.052 0.019
## Cumulative Var 0.374 0.426 0.445

print(modelo_afe$communality)

##      phq1      phq2      phq3      phq4      phq5      phq6      phq7      phq8 
## 0.4150915 0.7015517 0.4485456 0.5597728 0.3570552 0.5242933 0.3703462 0.3479024 
##      phq9 
## 0.2820829

La solución de cuatro factores presenta un buen ajuste y una estructura interpretable (RMSR = 0.025; índice de ajuste = 0.994). PA1 y PA2 muestran mayor precisión, PA3 una precisión aceptable y PA4 menor estabilidad, por lo que debe interpretarse con cautela.

fa.diagram(modelo_afe)

La rotación factorial facilita la interpretación al simplificar las cargas y puede ser ortogonal (Varimax) u oblicua (Oblimin, Quartimin, Promax), según se asuma independencia o correlación entre factores. En rotaciones oblicuas deben revisarse las matrices de patrones y estructura y las correlaciones entre factores, considerando valores > 0.80 como posible redundancia.

modelo_afe <- fa(matriz_cor, nfactors = 3, fm = "pa", rotate = "varimax")
salida <- capture.output(print(modelo_afe, cut = 0, digits = 3))
salida

##  [1] "Factor Analysis using method =  pa"                                       
##  [2] "Call: fa(r = matriz_cor, nfactors = 3, rotate = \"varimax\", fm = \"pa\")"
##  [3] "Standardized loadings (pattern matrix) based upon correlation matrix"     
##  [4] "       PA2   PA3   PA1    h2    u2  com"                                  
##  [5] "phq1 0.395 0.246 0.445 0.415 0.585 2.56"                                  
##  [6] "phq2 0.334 0.473 0.605 0.702 0.298 2.50"                                  
##  [7] "phq3 0.625 0.199 0.135 0.449 0.551 1.30"                                  
##  [8] "phq4 0.694 0.160 0.230 0.560 0.440 1.33"                                  
##  [9] "phq5 0.503 0.264 0.186 0.357 0.643 1.82"                                  
## [10] "phq6 0.227 0.560 0.400 0.524 0.476 2.18"                                  
## [11] "phq7 0.362 0.450 0.193 0.370 0.630 2.31"                                  
## [12] "phq8 0.290 0.508 0.074 0.348 0.652 1.64"                                  
## [13] "phq9 0.083 0.482 0.207 0.282 0.718 1.42"                                  
## [14] ""                                                                         
## [15] "                        PA2   PA3   PA1"                                  
## [16] "SS loadings           1.666 1.425 0.916"                                  
## [17] "Proportion Var        0.185 0.158 0.102"                                  
## [18] "Cumulative Var        0.185 0.343 0.445"                                  
## [19] "Proportion Explained  0.416 0.356 0.229"                                  
## [20] "Cumulative Proportion 0.416 0.771 1.000"                                  
## [21] ""                                                                         
## [22] "Mean item complexity =  1.9"                                              
## [23] "Test of the hypothesis that 3 factors are sufficient."                    
## [24] ""                                                                         
## [25] "df null model =  36  with the objective function =  2.551"                
## [26] "df of  the model are 12  and the objective function was  0.02 "           
## [27] ""                                                                         
## [28] "The root mean square of the residuals (RMSR) is  0.014 "                  
## [29] "The df corrected root mean square of the residuals is  0.024 "            
## [30] ""                                                                         
## [31] "Fit based upon off diagonal values = 0.999"                               
## [32] "Measures of factor score adequacy             "                           
## [33] "                                                    PA2   PA3    PA1"     
## [34] "Correlation of (regression) scores with factors   0.797 0.726  0.671"     
## [35] "Multiple R square of scores with factors          0.636 0.526  0.450"     
## [36] "Minimum correlation of possible factor scores     0.272 0.053 -0.101"

print(modelo_afe$loadings, cutoff = 0.30)

## 
## Loadings:
##      PA2   PA3   PA1  
## phq1 0.395       0.445
## phq2 0.334 0.473 0.605
## phq3 0.625            
## phq4 0.694            
## phq5 0.503            
## phq6       0.560 0.400
## phq7 0.362 0.450      
## phq8       0.508      
## phq9       0.482      
## 
##                  PA2   PA3   PA1
## SS loadings    1.666 1.425 0.916
## Proportion Var 0.185 0.158 0.102
## Cumulative Var 0.185 0.343 0.445

print(modelo_afe$communality)

##      phq1      phq2      phq3      phq4      phq5      phq6      phq7      phq8 
## 0.4150915 0.7015517 0.4485456 0.5597728 0.3570552 0.5242933 0.3703462 0.3479024 
##      phq9 
## 0.2820829

fa.diagram(modelo_afe)

Paso 7: Refinamiento del instrumento

Se recomienda depurar los ítems con cargas cruzadas, cargas débiles o problemas de codificación. Finalmente, se debe repetir el EFA con los ítems depurados y reevaluar KMO, Bartlett, ajuste y confiabilidad mediante alfa de Cronbach.

Es una técnica estadística multivariada y exploratoria que permite agrupar individuos u objetos en grupos homogéneos. Su objetivo es formar grupos cuyos elementos sean lo más similares posible dentro de cada grupo y lo más diferentes posible entre grupos.

A diferencia de la clasificación supervisada, los grupos no están previamente definidos, sino que se identifican a partir de las características observadas.

Los resultados deben interpretarse como estructuras descriptivas o exploratorias y no como relaciones causales.

El análisis puede organizarse en nueve etapas principales:

1. Definición del objetivo y unidad de análisis.
2. Selección de variables.
3. Análisis exploratorio de los datos.
4. Selección de la medida de distancia o disimilitud.
5. Selección del algoritmo de agrupamiento.
6. Determinación del número de conglomerados.
7. Evaluación de calidad y estabilidad.
8. Caracterización de los conglomerados.
9. Interpretación de la solución.

Antes del agrupamiento se debe definir:

• ¿Qué se desea identificar?
• ¿Qué elementos serán agrupados?

La unidad de análisis puede corresponder a:

• Individuos.
• Pacientes.
• Hogares.
• Instituciones.
• Regiones.
• Otros objetos de interés.

Las variables deben ser pertinentes al objetivo del estudio y aportar información para diferenciar las observaciones.

Su selección debe considerar:

• Naturaleza y escala de medición.
• Variabilidad.
• Posible redundancia entre variables.
• Presencia de valores atípicos.
• Cantidad y patrón de datos faltantes.

Se dice que es descriptivo y exploratorio, ya que las variables incluidas determinan directamente la estructura de los conglomerados.

Antes de calcular las distancias, es necesario realizar un análisis exploratorio que incluya:

• Distribución de las variables.
• Valores atípicos.
• Datos faltantes.
• Necesidad de transformación o estandarización.
• Relaciones entre variables cuantitativas.

Cuando las variables presentan escalas diferentes, se recomienda estandarizarlas para evitar que aquellas con mayor magnitud dominen el cálculo de las distancias.

Las distancias y disimilitudes cuantifican qué tan diferentes son dos observaciones.

Su elección depende de:

• Naturaleza de las variables.
• Escala de medición.
• Estructura esperada de los datos.
• Algoritmo de agrupamiento.

Para datos mixtos, una alternativa especialmente útil es la disimilitud de Gower.

Los métodos pueden clasificarse en tres grandes enfoques:

Jerárquicos

Construyen una estructura jerárquica mediante fusiones o divisiones sucesivas.

Particionamiento

Dividen directamente las observaciones en un número determinado de grupos (k-means)

Basados en densidad

Identifican regiones de alta densidad y pueden clasificar observaciones como ruido (DBSCAN)

El agrupamiento jerárquico construye una estructura de relaciones entre las observaciones.

Los resultados suelen representarse mediante un Dendrograma

Existen dos estrategias principales:

• Aglomerativa: de abajo hacia arriba.
• Divisiva: de arriba hacia abajo.

El procedimiento es:

1. Cada observación comienza como un conglomerado.
2. Se calculan las distancias entre conglomerados.
3. Se fusionan los dos más próximos.
4. Se actualizan las distancias.
5. Se repite el proceso hasta obtener un único conglomerado.

Es el enfoque jerárquico más utilizado.

El procedimiento sigue una estrategia de arriba hacia abajo:

1. Todas las observaciones comienzan en un único conglomerado.
2. El conjunto se divide progresivamente.
3. Las divisiones continúan hasta alcanzar el nivel definido por el algoritmo.

El método de enlace define cómo se resume la proximidad entre dos conglomerados.

• Enlace simple: utiliza la menor distancia. Puede generar grupos alargados por efecto de encadenamiento.
• Enlace completo: utiliza la mayor distancia. Favorece grupos compactos, aunque puede ser sensible a valores atípicos.
• Enlace promedio: considera todas las distancias y representa un equilibrio entre ambos.
• Enlace de centroide: utiliza la distancia entre los centros de los conglomerados y es principalmente apropiado para variables cuantitativas.
• Ward: busca minimizar el aumento de la variabilidad interna, favoreciendo conglomerados compactos y homogéneos, aunque puede ser sensible a valores atípicos.

Los métodos de particionamiento asignan directamente las observaciones a grupos sin construir un dendrograma y, generalmente, requieren definir previamente el número de conglomerados.

K-means es uno de los métodos más utilizados:

1. Asigna cada observación al centroide más cercano.
2. Recalcula los centroides.
3. Reasigna las observaciones.
4. Repite el proceso hasta alcanzar la convergencia.

Como la solución puede depender de los centroides iniciales, se recomienda utilizar múltiples inicializaciones (nstart).

Los métodos basados en densidad forman conglomerados a partir de regiones con alta concentración de observaciones, separadas por zonas de menor densidad.

A diferencia de K-means, pueden:

• Identificar grupos de formas irregulares.
• Clasificar algunas observaciones como ruido.
• No requerir necesariamente la especificación previa del número de grupos.

DBSCAN utiliza dos parámetros principales:

• ε: radio de vecindad.
• MinPts: número mínimo de observaciones en la vecindad.

Clasifica las observaciones como:

• Puntos núcleo.
• Puntos frontera.
• Ruido.

Entre sus ventajas se encuentran que no requiere definir previamente el número de conglomerados, permite detectar formas irregulares e identifica explícitamente el ruido.

Su desempeño depende de la elección de ε y MinPts, y puede presentar dificultades cuando los grupos tienen densidades diferentes.

Las medidas de similitud cuantifican qué tan parecidas son dos observaciones y cumplen una función diferente de las distancias o disimilitudes utilizadas para formar conglomerados.

• Jaccard: útil para variables binarias asimétricas; valores entre 0 y 1.
• Dice: similar a Jaccard, pero otorga mayor peso a las coincidencias positivas.
• Phi: mide la asociación entre dos variables binarias; valores entre −1 y 1.
• V de Cramer: evalúa la intensidad de asociación entre variables categóricas; valores entre 0 y 1.
• Coeficiente de contingencia: mide asociación entre variables categóricas, aunque su máximo depende del tamaño de la tabla.

Evalúa la calidad, cohesión, separación y estabilidad de los conglomerados. Puede ser interna, basada en los datos, o externa, comparando con una clasificación independiente.

• Silhouette: valores cercanos a 1 indican buena asignación; valores negativos, posible mala clasificación.
• Dunn: valores altos indican grupos compactos y bien separados.
• Davies-Bouldin: valores bajos indican mejor agrupamiento.
• Calinski-Harabasz: valores altos indican mayor separación y menor variabilidad interna.

Cuando existe una clasificación independiente, pueden utilizarse medidas como:

• Jaccard
• Rand
• Rand ajustado
• Información mutua
• Información mutua ajustada
• V de Cramer

Estas medidas evalúan la concordancia entre la solución obtenida y una clasificación externa. La clasificación de referencia debe considerarse independiente y no necesariamente como una verdad absoluta.

La elección del número de conglomerados (k) es fundamental, ya que influye en la estructura, interpretación y utilidad de la solución.

No existe un criterio único que sea óptimo en todos los casos.

Se recomienda combinar:

• Medidas de validación.
• Estabilidad.
• Interpretabilidad.
• Conocimiento sustantivo del problema.

En K-means, evalúa la variabilidad interna para distintos valores de k. Al aumentar k, la variabilidad disminuye. Se selecciona el punto donde la mejora se reduce notablemente, denominado “codo”. Su identificación puede ser subjetiva.

• Silhouette: valores mayores indican mejor cohesión y separación.
• Davies-Bouldin: valores menores indican una mejor estructura.
• Calinski-Harabasz: valores mayores indican mayor separación y menor dispersión interna.
• CCC: valores altos o picos locales pueden respaldar determinadas soluciones, aunque debe interpretarse junto con otros criterios.

La selección de k debe considerar un rango de valores, comparando diferentes índices de validación y evaluando la estabilidad de las soluciones, el tamaño de los conglomerados y la interpretabilidad de los grupos.

Una vez seleccionada la solución, los conglomerados se caracterizan según las variables utilizadas, identificando las características que los diferencian, los perfiles predominantes y posibles grupos pequeños. La interpretación debe tener sentido sustantivo, recordando que el análisis de conglomerados identifica perfiles o estructuras, no relaciones causales.

En este ejemplo se utilizará K-means para identificar perfiles de participantes a partir de características demográficas, antropométricas, cardiovasculares, metabólicas y de salud mental.

Posteriormente se comparará la distribución de un problema de sueño entre los conglomerados.

Objetivo

Identificar perfiles homogéneos de participantes y explorar si el problema de sueño se distribuye de manera diferente entre los perfiles.

Para construir los conglomerados utilizaremos: Edad, IMC, Circunferencia de cintura, Presión arterial Sistólica, HbA1c, Colesterol, HDL, Puntaje PHQ.

Nota problema_sueno_consulta no se utilizará para construir los conglomerados, sino como una variable suplementaria.

Paso 2. Preparar los datos

set.seed(57123)

# Seleccionar variables para construir los clústeres
datos_cluster <- base %>%
  dplyr::select(edad,imc,cintura,pas1,pad1,hba1c,colesterol,hdl,phq_total) 

Paso 3. Estandarizar las variables

Las variables están expresadas en escalas diferentes, por lo tanto, se estandarizan antes de aplicar K-means.

datos_scaled <- scale(datos_cluster)

No debemos asumir un número específico de grupos. Se recomienda evaluar diferentes valores de k.

Método del codo

fviz_nbclust(datos_scaled, kmeans, method = "wss") +
  labs(title = "Método del codo", x = "Número de clústeres (k)",
    y = "Suma de cuadrados intra-clúster" )

El gráfico permite observar cuánto disminuye la variabilidad interna al aumentar el número de grupos.

También podemos utilizar el índice de Silhouette:

fviz_nbclust(datos_scaled, kmeans, method = "silhouette") +
  labs( title = "Método de la silueta", x = "Número de clústeres (k)",
    y = "Ancho promedio de la silueta")

El valor de k sugerido por un índice no debe aceptarse automáticamente. Se debe contrastar con otros criterios y con la interpretabilidad de los grupos.

Para el ejemplo se comienza con k = 3.

k <- 3
modelo_kmeans <- kmeans(datos_scaled, centers = k, nstart = 25)
modelo_kmeans

## K-means clustering with 3 clusters of sizes 1216, 1060, 1234
## 
## Cluster means:
##          edad        imc    cintura       pas1       pad1      hba1c
## 1 -0.84593241 -0.5432874 -0.6992812 -0.7390304 -0.2815286 -0.4693647
## 2  0.02946535  1.0916659  1.1163534  0.1149074  0.1632981  0.4116067
## 3  0.80828245 -0.4023730 -0.2698612  0.6295455  0.1371497  0.1089500
##    colesterol         hdl   phq_total
## 1 -0.24872408  0.08542631 -0.07173102
## 2 -0.09858894 -0.55894850  0.22018753
## 3  0.32978344  0.39595382 -0.11845532
## 
## Clustering vector:
##    [1] 1 1 3 3 1 2 3 3 3 2 1 1 1 3 3 1 3 1 3 3 3 1 3 2 3 1 2 3 1 1 2 1 1 3 2 1 2
##   [38] 3 3 2 3 1 2 1 1 2 2 1 3 3 3 3 3 2 2 1 3 1 1 1 2 3 1 3 3 1 2 2 3 1 3 3 3 1
##   [75] 1 1 2 3 2 2 2 2 1 1 1 2 1 2 1 3 2 3 1 2 3 2 1 3 1 3 1 3 3 2 2 1 3 2 3 1 1
##  [112] 3 1 2 1 2 3 1 3 1 2 1 1 3 1 1 3 1 3 3 1 1 3 3 3 1 1 1 3 3 3 2 3 3 1 3 1 2
##  [149] 1 2 3 2 3 3 2 2 3 2 3 2 2 1 3 3 3 2 3 2 1 3 1 3 1 3 3 3 1 2 2 1 3 2 3 1 3
##  [186] 2 3 3 1 3 1 1 3 1 3 3 1 3 2 2 3 1 1 2 2 3 1 2 3 3 2 3 1 3 1 3 2 2 1 1 3 1
##  [223] 2 1 2 1 2 2 3 3 1 3 1 2 3 3 1 2 2 3 1 1 1 3 1 3 2 3 3 1 1 1 1 1 2 3 3 2 2
##  [260] 3 3 2 2 2 1 2 2 3 3 2 1 2 1 3 2 1 2 1 2 1 2 1 3 2 3 3 1 3 3 1 2 1 2 2 3 3
##  [297] 3 2 2 2 2 2 2 1 1 1 2 2 2 1 1 1 2 1 1 3 3 1 2 2 3 1 3 1 3 2 1 1 1 3 1 1 3
##  [334] 3 2 1 1 1 3 3 3 3 3 3 3 2 2 2 2 3 1 2 3 2 1 3 3 3 1 3 1 1 2 3 1 3 2 2 2 1
##  [371] 1 1 2 1 1 3 3 3 2 1 2 3 2 3 3 1 1 3 1 1 2 2 3 1 3 3 1 2 1 1 2 3 1 1 3 3 3
##  [408] 2 1 2 3 1 1 2 1 1 3 1 2 3 2 3 3 2 1 2 2 1 1 2 1 1 1 3 2 3 1 1 1 1 3 1 3 1
##  [445] 1 1 3 1 3 3 3 2 2 1 1 2 2 2 2 2 3 2 2 2 3 2 1 1 3 3 2 2 3 1 1 2 2 2 1 3 2
##  [482] 1 3 3 1 3 2 1 1 3 1 1 2 3 3 3 1 3 1 2 1 3 3 1 3 1 2 2 1 3 3 3 1 3 3 2 2 3
##  [519] 2 1 3 1 2 3 2 1 2 1 2 3 3 1 3 2 3 2 3 3 1 3 2 3 1 2 1 2 2 3 3 2 3 3 1 3 3
##  [556] 1 2 1 3 3 1 1 1 3 2 3 3 2 3 3 3 2 3 2 1 2 1 2 3 1 2 2 2 3 1 3 1 1 1 3 2 1
##  [593] 3 1 3 1 1 1 2 3 1 3 1 1 2 2 2 2 1 3 1 2 1 1 3 3 3 2 1 1 2 2 2 3 1 1 1 3 1
##  [630] 3 2 2 3 3 3 2 1 2 1 3 2 1 1 2 2 2 3 3 2 2 3 1 3 1 3 1 1 3 3 1 1 3 2 3 3 1
##  [667] 1 2 1 3 1 1 2 1 3 2 1 1 1 3 2 3 2 3 2 3 3 3 3 3 1 1 1 1 2 1 3 3 3 2 1 3 2
##  [704] 2 2 1 1 3 3 3 1 3 1 3 2 2 3 3 2 3 2 3 1 2 3 1 1 2 1 1 2 1 1 1 3 3 3 1 2 3
##  [741] 3 3 3 2 3 1 2 2 1 1 1 1 1 3 1 3 2 2 2 1 3 1 1 3 3 2 1 2 2 1 1 1 2 1 2 3 3
##  [778] 3 2 3 1 2 1 2 1 2 2 3 3 1 2 2 1 2 3 1 3 3 2 1 3 3 3 3 1 2 1 1 2 3 2 2 2 3
##  [815] 2 3 3 3 1 1 2 1 1 3 2 1 3 2 2 3 1 3 2 3 1 3 3 2 3 2 2 1 2 3 3 3 2 1 2 1 3
##  [852] 2 3 2 3 1 3 2 2 3 1 2 1 3 1 2 1 2 3 1 1 3 2 2 2 2 3 2 3 1 1 3 2 3 2 1 2 2
##  [889] 1 1 3 3 3 3 1 1 3 3 1 3 3 1 3 3 3 2 3 3 1 1 1 2 3 2 3 2 3 3 1 3 1 1 1 2 1
##  [926] 2 2 2 2 1 2 3 2 1 1 2 3 2 1 2 2 3 3 1 3 1 3 2 3 3 2 1 1 3 3 2 1 2 2 2 1 2
##  [963] 1 3 3 3 3 2 3 3 2 3 3 2 3 2 1 1 3 2 1 2 2 2 1 1 1 1 2 2 3 3 1 1 2 2 2 1 3
## [1000] 3 1 1 3 1 1 2 2 2 3 2 2 3 1 3 1 1 3 2 3 2 3 1 1 2 3 2 1 1 3 1 1 3 1 1 2 1
## [1037] 3 2 2 2 3 3 1 3 1 2 2 3 3 2 3 1 2 1 1 3 1 3 1 1 2 1 3 2 3 2 1 1 1 3 3 1 2
## [1074] 1 2 2 2 3 1 2 1 1 2 3 2 3 1 1 3 3 2 2 3 3 1 1 3 2 1 2 3 1 1 3 1 3 3 1 3 3
## [1111] 2 3 3 3 1 2 1 3 2 3 1 3 1 3 3 1 1 2 2 1 3 3 2 2 2 1 1 3 2 2 3 2 2 1 2 3 1
## [1148] 1 2 1 3 1 1 1 1 2 2 3 1 3 3 1 3 1 2 2 3 3 1 1 1 2 3 1 1 2 1 1 1 2 3 1 2 2
## [1185] 3 3 1 3 1 1 3 3 2 2 1 1 2 2 1 2 2 2 2 2 3 1 1 2 3 3 3 1 3 2 3 3 2 2 2 3 2
## [1222] 2 3 1 3 1 2 1 3 2 3 2 3 3 3 3 3 3 1 2 1 1 3 3 2 3 3 3 1 1 3 3 1 1 3 2 2 3
## [1259] 2 1 2 2 2 3 1 2 3 3 1 2 1 1 2 1 3 2 2 3 2 1 2 3 3 1 1 1 1 1 3 3 1 2 2 1 2
## [1296] 2 1 3 1 1 3 3 3 1 1 2 1 1 1 1 1 3 3 1 1 3 1 1 3 3 2 1 1 1 2 2 1 3 2 1 3 3
## [1333] 3 2 3 1 3 1 3 1 3 2 2 2 1 1 1 3 1 1 2 2 3 3 2 3 3 3 1 1 3 1 3 2 1 2 2 1 3
## [1370] 3 1 3 3 1 3 2 3 2 1 3 3 1 3 2 1 1 1 3 1 1 2 1 3 2 1 1 3 1 2 1 2 3 1 1 2 2
## [1407] 1 2 2 1 3 2 1 3 3 2 3 1 1 3 3 2 3 3 3 1 3 1 2 2 1 1 1 2 3 3 2 2 2 1 2 2 3
## [1444] 2 3 3 3 1 3 2 2 3 3 3 1 2 2 1 2 1 3 3 3 3 2 1 1 2 2 1 1 1 3 1 1 2 1 3 1 2
## [1481] 2 3 3 1 2 2 3 3 3 2 1 3 1 1 1 3 3 2 1 1 3 1 3 1 2 1 3 1 3 3 1 1 2 2 1 1 3
## [1518] 1 1 1 3 2 1 3 1 1 3 3 3 3 1 3 2 3 1 2 3 1 2 2 2 3 2 3 2 2 1 1 3 2 2 2 2 2
## [1555] 3 3 2 1 1 2 3 1 3 2 2 1 3 2 2 3 1 1 1 3 3 1 3 2 2 2 3 2 2 3 3 2 2 3 1 1 3
## [1592] 1 3 2 1 3 3 1 1 2 3 1 3 2 3 2 3 3 3 3 3 1 2 1 2 3 2 1 1 3 3 3 2 1 3 1 1 1
## [1629] 1 1 1 3 3 3 1 3 2 1 1 2 2 1 2 3 2 2 3 3 3 3 1 1 2 1 1 2 1 3 1 2 3 1 1 2 3
## [1666] 1 2 3 2 3 1 1 1 3 1 2 2 2 1 1 3 2 2 3 1 3 3 2 3 3 1 1 3 3 3 2 2 3 1 3 2 1
## [1703] 1 1 2 3 3 1 3 3 1 3 1 1 2 2 3 1 2 2 3 3 1 3 3 2 3 2 3 2 3 2 1 2 1 2 1 2 2
## [1740] 1 3 1 2 2 1 3 3 2 3 1 2 3 3 2 1 3 1 2 1 3 3 2 2 3 3 2 1 2 1 3 1 1 3 2 3 1
## [1777] 3 2 2 1 1 1 2 1 2 2 3 2 1 1 3 3 1 3 1 1 1 1 3 1 1 3 1 2 1 2 2 2 3 2 3 1 2
## [1814] 3 2 1 3 3 2 3 3 1 1 3 2 3 3 3 3 3 3 1 3 3 2 1 3 2 3 2 3 3 3 2 3 3 2 3 2 1
## [1851] 3 3 3 1 3 2 1 3 3 1 3 2 2 1 1 3 3 2 2 3 3 1 3 2 3 1 1 3 2 1 3 3 1 3 2 3 2
## [1888] 3 1 3 2 3 1 2 1 2 1 2 3 3 2 1 1 2 1 2 3 1 2 3 1 2 1 1 2 1 2 2 3 2 1 2 2 1
## [1925] 1 3 2 2 2 1 2 1 3 1 2 1 1 3 2 2 1 3 3 3 2 3 3 3 2 2 3 1 3 3 1 1 2 1 1 3 3
## [1962] 2 1 2 1 2 1 3 1 1 1 1 3 3 1 3 2 3 2 2 3 2 1 2 1 3 2 2 2 3 2 2 1 3 1 1 1 2
## [1999] 2 2 1 3 2 1 2 1 1 3 3 3 3 2 1 1 3 1 1 1 3 3 3 2 1 3 2 3 3 1 1 2 2 3 3 1 3
## [2036] 1 3 1 3 1 1 3 2 2 3 1 1 3 1 2 3 1 1 3 3 1 1 3 1 3 3 1 3 3 1 2 3 3 2 1 1 3
## [2073] 1 1 2 1 3 1 2 1 1 1 2 3 1 1 1 1 2 2 2 3 1 2 3 1 1 2 2 3 1 2 2 3 1 3 1 2 1
## [2110] 2 1 2 3 1 1 1 3 1 3 1 1 3 3 3 1 3 2 2 2 1 1 1 2 2 3 3 2 3 1 3 3 1 1 2 2 2
## [2147] 3 3 3 1 2 3 2 2 2 1 2 3 3 3 3 1 2 2 1 1 2 3 1 2 1 2 1 1 1 1 2 2 2 3 3 2 1
## [2184] 3 2 3 1 3 3 1 1 3 1 1 3 2 2 1 2 3 2 2 2 3 2 3 3 3 2 2 2 2 2 3 1 3 1 1 2 3
## [2221] 1 3 1 1 2 1 1 1 3 2 3 2 1 2 2 1 1 2 2 1 2 1 1 3 3 3 1 3 2 3 3 2 1 3 3 3 1
## [2258] 1 2 2 2 3 3 3 1 3 2 2 1 3 1 3 3 3 3 3 3 3 3 1 3 2 1 3 2 2 3 3 3 2 2 1 3 2
## [2295] 2 2 1 1 2 2 3 1 3 2 3 1 1 2 3 3 3 3 1 2 2 1 2 2 3 2 2 1 2 2 1 3 1 2 2 3 1
## [2332] 1 3 3 3 1 1 1 2 3 2 3 3 1 1 2 2 1 3 1 1 1 2 3 3 1 1 1 1 1 1 3 2 3 1 3 1 2
## [2369] 3 3 3 1 3 2 2 3 1 1 2 1 2 2 1 1 2 3 1 1 1 2 1 2 1 1 1 3 2 2 1 2 2 3 3 2 3
## [2406] 2 1 3 1 3 3 2 3 1 1 2 1 2 2 1 1 2 2 1 1 2 1 1 1 3 2 3 2 2 2 2 2 1 1 3 3 2
## [2443] 3 2 1 3 3 1 1 1 2 3 2 3 2 1 2 1 1 2 1 1 1 1 3 3 2 1 3 1 2 2 1 2 3 1 2 2 1
## [2480] 3 3 1 2 3 3 1 3 1 2 1 2 1 2 2 2 1 3 2 3 2 3 1 2 2 1 1 1 1 1 3 3 1 1 3 2 3
## [2517] 3 1 3 3 3 3 1 2 2 2 3 3 2 2 2 3 2 2 3 1 2 1 1 2 1 1 1 2 3 2 1 2 2 1 3 3 1
## [2554] 2 3 3 2 3 1 3 2 3 3 3 1 3 2 3 1 1 2 2 3 1 2 2 1 3 3 2 2 2 2 3 2 3 2 1 3 3
## [2591] 3 2 1 2 1 3 1 1 2 1 1 3 2 2 1 3 2 1 1 1 2 1 1 1 1 3 3 2 1 1 2 2 1 1 3 1 1
## [2628] 2 3 2 1 1 1 3 2 2 2 2 1 3 1 3 3 1 3 1 1 3 3 3 3 3 3 1 3 3 3 3 3 3 2 2 3 2
## [2665] 3 1 3 2 3 3 3 2 3 1 2 2 3 3 1 3 2 2 3 2 2 3 3 3 1 3 2 1 3 1 3 2 3 2 2 1 1
## [2702] 2 3 3 1 3 1 3 2 1 3 3 3 3 2 1 2 3 1 2 3 2 1 1 3 3 3 2 1 2 2 1 3 1 1 1 1 2
## [2739] 2 1 1 2 3 1 1 1 3 3 1 3 1 3 2 2 2 2 3 3 1 1 1 1 1 3 1 1 2 3 1 2 1 3 3 2 1
## [2776] 3 2 1 1 3 2 2 1 1 2 2 3 2 1 2 2 1 2 2 2 2 2 3 3 1 1 3 1 3 3 1 3 1 3 2 3 1
## [2813] 1 3 1 3 3 3 1 2 2 1 2 3 3 1 1 3 1 3 3 2 1 3 1 2 3 1 2 3 2 3 3 1 1 2 2 2 1
## [2850] 2 3 2 3 1 1 2 3 1 2 3 2 3 3 1 3 3 3 3 2 3 2 2 3 3 1 1 3 3 2 3 1 1 1 1 3 1
## [2887] 3 3 2 3 1 3 2 2 1 1 3 3 1 1 1 2 1 1 2 1 2 2 2 1 1 2 2 3 2 2 3 2 2 2 3 1 3
## [2924] 3 1 3 3 3 3 2 2 1 2 1 2 1 1 3 1 1 3 1 2 3 1 2 1 1 2 2 1 3 1 3 3 3 1 2 1 1
## [2961] 1 1 1 1 3 1 3 1 2 1 3 2 3 1 3 3 2 1 3 3 2 2 2 2 1 2 1 1 3 3 3 2 1 2 1 3 3
## [2998] 1 3 1 2 3 1 3 1 2 3 3 3 1 3 2 1 3 2 2 2 2 1 3 3 3 3 2 2 3 2 1 1 1 2 1 2 1
## [3035] 2 3 2 3 3 1 1 3 3 1 1 3 2 2 2 2 1 1 3 1 3 2 3 2 2 3 2 1 2 2 3 3 3 1 3 2 3
## [3072] 1 3 3 1 1 3 2 1 3 3 1 3 1 1 1 3 3 2 1 3 2 3 3 1 1 3 1 2 3 1 2 2 2 3 2 1 3
## [3109] 2 1 1 3 3 3 2 3 3 2 2 2 1 3 3 2 2 1 2 2 3 2 1 1 1 2 3 3 3 3 2 3 2 3 1 1 3
## [3146] 2 2 1 3 3 2 2 2 3 1 1 2 3 1 2 2 1 1 1 2 3 3 2 1 1 1 3 1 1 3 2 3 3 3 3 3 3
## [3183] 2 2 3 3 1 2 1 2 2 1 2 1 2 2 3 1 3 3 2 1 3 1 3 1 1 1 1 2 2 1 2 1 1 1 1 3 1
## [3220] 3 2 1 3 2 3 1 2 3 1 1 3 3 2 2 2 1 2 1 2 3 1 3 3 3 1 1 2 2 1 1 2 3 1 1 3 2
## [3257] 1 3 3 2 3 3 2 2 3 1 1 2 1 1 1 2 2 3 1 1 3 1 1 3 2 2 3 2 2 1 3 2 3 2 2 2 2
## [3294] 1 3 1 3 3 2 1 3 1 1 1 1 1 3 1 3 1 3 2 1 2 2 2 1 2 1 3 3 1 2 3 1 3 2 3 2 3
## [3331] 3 2 1 2 3 2 3 2 3 1 1 2 1 2 3 2 1 2 3 1 3 3 1 3 3 2 1 3 3 3 2 1 3 1 2 3 2
## [3368] 1 1 3 1 3 3 1 2 1 3 1 2 3 2 1 2 3 1 3 2 1 2 1 1 1 1 1 1 3 3 2 2 3 1 3 1 1
## [3405] 3 1 2 2 1 3 3 1 3 1 3 3 2 2 3 2 2 3 2 1 3 3 3 2 1 3 1 2 3 1 2 2 3 2 1 2 2
## [3442] 1 1 3 3 1 3 2 2 2 1 3 2 3 1 1 1 1 3 1 3 1 1 1 3 2 2 1 3 1 1 2 3 3 3 2 2 2
## [3479] 3 2 1 3 3 2 2 3 1 1 3 2 3 1 1 2 3 2 3 2 3 1 1 3 2 1 2 3 1 3 1 2
## 
## Within cluster sum of squares by cluster:
## [1] 5670.136 8626.397 9174.344
##  (between_SS / total_SS =  25.7 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Se utiliza nstart = 25 para ejecutar el algoritmo desde múltiples configuraciones iniciales y reducir la dependencia de una única inicialización.

modelo_kmeans$size

## [1] 1216 1060 1234

Este resultado permite comprobar cuántas observaciones pertenecen a cada grupo. Un conglomerado extremadamente pequeño puede ser una señal para revisar la solución.

modelo_kmeans$centers

##          edad        imc    cintura       pas1       pad1      hba1c
## 1 -0.84593241 -0.5432874 -0.6992812 -0.7390304 -0.2815286 -0.4693647
## 2  0.02946535  1.0916659  1.1163534  0.1149074  0.1632981  0.4116067
## 3  0.80828245 -0.4023730 -0.2698612  0.6295455  0.1371497  0.1089500
##    colesterol         hdl   phq_total
## 1 -0.24872408  0.08542631 -0.07173102
## 2 -0.09858894 -0.55894850  0.22018753
## 3  0.32978344  0.39595382 -0.11845532

Los centros están expresados en unidades estandarizadas.

• Valores positivos → por encima del promedio.
• Valores negativos → por debajo del promedio.
• Valores cercanos a 0 → próximos al promedio.

Esto permite construir un perfil de cada conglomerado.

Primero calculamos las distancias entre las observaciones:

distancias <- dist(datos_scaled)

Luego calculamos la silueta:

sil <- silhouette(modelo_kmeans$cluster, distancias)
summary(sil)

## Silhouette of 3510 units in 3 clusters from silhouette.default(x = modelo_kmeans$cluster, dist = distancias) :
##  Cluster sizes and average silhouette widths:
##       1216       1060       1234 
## 0.26809744 0.09150828 0.08400072 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.13680  0.05223  0.14148  0.15005  0.24090  0.45158

fviz_silhouette(sil) +
  labs(
    title = "Calidad de los clústeres mediante silueta"
  )

##   cluster size ave.sil.width
## 1       1 1216          0.27
## 2       2 1060          0.09
## 3       3 1234          0.08

Una solución adecuada debería presentar:

• Valores promedio de Silhouette relativamente altos.
• Pocas observaciones con valores negativos.
• Grupos con buena separación.

fviz_cluster(modelo_kmeans, data = datos_scaled,
  geom = "point", ellipse.type = "convex",
  main = "Perfiles identificados mediante K-means")

Este gráfico proporciona una representación bidimensional de la estructura de los conglomerados.

Asignamos a cada observación el conglomerado correspondiente:

datos_resultado <- datos_cluster %>%
  mutate(cluster = factor(modelo_kmeans$cluster))

Podemos resumir las características de cada grupo:

resumen_cluster <- datos_resultado %>%
  group_by(cluster) %>%
  summarise(n = n(), edad = median(edad, na.rm = TRUE),
    imc = median(imc, na.rm = TRUE),
    cintura = median(cintura, na.rm = TRUE),
    pas1 = median(pas1, na.rm = TRUE),
    pad1 = median(pad1, na.rm = TRUE),
    hba1c = median(hba1c, na.rm = TRUE),
    colesterol = median(colesterol, na.rm = TRUE),
    hdl = median(hdl, na.rm = TRUE),
    phq_total = median(phq_total, na.rm = TRUE))
resumen_cluster

## # A tibble: 3 × 11
##   cluster     n  edad   imc cintura  pas1  pad1 hba1c colesterol   hdl phq_total
##   <fct>   <int> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>      <dbl> <dbl>     <dbl>
## 1 1        1216    34  25.7    89.6   112    70   5.3        175    53         2
## 2 2        1060    52  36.1   118.    126    74   5.8        183    43         3
## 3 3        1234    65  26.8    97     136    74   5.7        200    57         1

A partir de esta tabla podemos construir perfiles como:

• Clúster 1: participantes relativamente jóvenes, con menor IMC y menor riesgo metabólico.
• Clúster 2: participantes de mayor edad, con mayor adiposidad y presión arterial.
• Clúster 3: participantes con mayor puntaje de síntomas depresivos y características metabólicas intermedias.

Una vez construidos y evaluados los conglomerados, podemos estudiar una variable externa. En este caso:

base_cluster <- base %>%
  dplyr::select(edad, imc, cintura, pas1, pad1, hba1c, colesterol,
                hdl, phq_total, problema_sueno_consulta) %>%
  na.omit() %>%
  mutate(
    cluster = factor(modelo_kmeans$cluster)
  )

Construimos una tabla de contingencia:

tabla_sueno <- table(
  Cluster = base_cluster$cluster,
  Problema_sueno = base_cluster$problema_sueno_consulta
)

tabla_sueno

##        Problema_sueno
## Cluster  No  Sí
##       1 968 248
##       2 641 419
##       3 870 364

Para obtener porcentajes dentro de cada clúster:

prop.table(tabla_sueno, margin = 1)

##        Problema_sueno
## Cluster        No        Sí
##       1 0.7960526 0.2039474
##       2 0.6047170 0.3952830
##       3 0.7050243 0.2949757

¿Qué proporción de participantes con problema de sueño existe dentro de cada perfil?

ggplot(base_cluster, aes(x = cluster,fill = problema_sueno_consulta)
) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Problema de sueño según perfil identificado",
    x = "Clúster", y = "Proporción", fill = "Problema de sueño") +
  theme_minimal()

El gráfico muestra la proporción de participantes con y sin problema de sueño dentro de cada conglomerado.

Finalmente, podemos evaluar si existe evidencia de asociación entre el conglomerado y el problema de sueño mediante una prueba de chi-cuadrado:

chisq.test(tabla_sueno)

## 
##  Pearson's Chi-squared test
## 
## data:  tabla_sueno
## X-squared = 99.954, df = 2, p-value < 2.2e-16

Interpretación

La prueba permite evaluar si la distribución del problema de sueño es independiente del conglomerado. Sin embargo, una asociación estadística no implica causalidad.

Para evaluar el desempeño predictivo de los modelos supervisados se siguió como protocolo la separación de los datos, ajuste del modelo, selección del punto de corte y evaluación en una muestra independiente. El objetivo es estimar el rendimiento del modelo fuera de la muestra utilizada para su construcción y reducir el riesgo de sobreajuste.

                         [BASE DE DATOS COMPLETA]
                                  │
                 ┌────────────────┴────────────────┐
                 ▼                                 ▼
          [ENTRENAMIENTO 70%]               [VALIDACIÓN 30%]
                 │                                 │
          Ajuste del modelo                 Predicciones
                 │                                 │
       Selección de variables                Evaluación
       y/o hiperparámetros                   out-of-sample
                 │                                 │
       Selección del umbral                    ┌───┴────┐
          óptimo (Youden)                      ▼        ▼
                 │                           AUC   Matriz de
                 └──────────────────────────►      confusión

La muestra de entrenamiento (70%) se utiliza exclusivamente para construir y optimizar el modelo, incluyendo la selección de variables, hiperparámetros y punto de corte óptimo mediante el índice de Youden.

Posteriormente, el modelo final se aplica sin modificaciones a la muestra de validación (30%), que permanece independiente del proceso de ajuste. El desempeño predictivo se evalúa mediante el AUC, matriz de confusión y métricas de clasificación, permitiendo valorar la capacidad de generalización del modelo.

Método 5. Regresión logística múltiple

La regresión logística múltiple es una extensión de la regresión logística simple. Se basa en los mismos principios que esta, pero ampliando el número de predictores (\(x_i\)), los cuales pueden ser tanto continuos como categóricos.

\[\ln\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_i x_i\]

\[\text{logit}(Y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_i x_i\]

El valor de la probabilidad de \(Y\) se obtiene aplicando la función inversa del logit (función sigmoide):

\[p(Y) = \frac{e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_i x_i}}{1 + e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_i x_i}}\]

A la hora de evaluar la validez y calidad de un modelo de regresión logística múltiple, se analiza tanto el modelo en su conjunto como los predictores que lo forman:

• Comparación con el modelo nulo: Se considera que el modelo es útil si demuestra una mejora significativa respecto al modelo nulo (el modelo que no incluye ningún predictor).
• Pruebas de significación global: Existen tres tests estadísticos principales que cuantifican esta mejora mediante la comparación de las devianzas o residuos:
  1. Likelihood Ratio Test (Prueba de Razón de Verosimilitud)
  2. Score Test
  3. Wald Test

Tomando como variable de desenlace la presencia de problemas de sueño y como predictores las variables sociodemográficas, clínicas y de estilo de vida.

# Asignando categorías de referencia para las variables categóricas
base_modelos$grupo_edad <- factor(base_modelos$grupo_edad,  levels = c("18-30", "31-45","46-60", "61+"))
base_modelos$grupo_edad <- relevel(base_modelos$grupo_edad,ref = "18-30")
base_modelos$hba1ccat <- relevel(base_modelos$hba1ccat, ref = "Normal")
base_modelos$hta <- factor(base_modelos$hta, levels = c("No", "Sí"))
base_modelos$hta <- relevel(base_modelos$hta,  ref = "No")
base_modelos$raza <- relevel(base_modelos$raza, ref="Blanco")
base_modelos$sexo <- relevel(base_modelos$sexo, ref = "Hombre")
base_modelos$educacion <- relevel(base_modelos$educacion, ref="Univ/Posg")
base_modelos$quintil_pobreza <- factor(base_modelos$quintil_pobreza, levels= c("Muy bajo", "Bajo", "Medio", "Alto", "Muy alto"))
base_modelos$quintil_pobreza <- relevel(base_modelos$quintil_pobreza, ref = "Muy bajo")
base_modelos$imccat <- relevel(base_modelos$imccat, ref = "Normal")
base_modelos$actividad_vigorosa <- relevel(base_modelos$actividad_vigorosa, ref = "No")
base_modelos$fumo100 <- relevel(base_modelos$fumo100, ref = "No")
base_modelos$alcohol_12m <- relevel(base_modelos$alcohol_12m, ref = "No")
base_modelos$problema_sueno_consulta <- relevel(base_modelos$problema_sueno_consulta, ref = "No")

Un problema que tiene el PHQ-9 es que una de sus preguntas está relacionada con el sueño, por lo que para evitar la colinealidad, se propone crear un nuevo puntaje total del PHQ-9 sin incluir la pregunta relacionada con el sueño.

items_phq <- c("phq1", "phq2", "phq4", "phq5","phq6", "phq7", "phq8", "phq9")
base_modelos$phq_sin_sueno <- rowSums(sapply(base_modelos[items_phq],function(x) as.numeric(x) - 1),na.rm = FALSE)

modelo_log <- modelo_log <- glm(problema_sueno_consulta ~ grupo_edad + raza + educacion + estado_civil + quintil_pobreza + actividad_vigorosa + alcohol_12m + fumo100 + imc + phq_sin_sueno, 
                                data = base_modelos, family = binomial(link = "logit"))

summary(modelo_log)

## 
## Call:
## glm(formula = problema_sueno_consulta ~ grupo_edad + raza + educacion + 
##     estado_civil + quintil_pobreza + actividad_vigorosa + alcohol_12m + 
##     fumo100 + imc + phq_sin_sueno, family = binomial(link = "logit"), 
##     data = base_modelos)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -3.517072   0.323035 -10.888  < 2e-16 ***
## grupo_edad31-45           0.529013   0.153377   3.449 0.000562 ***
## grupo_edad46-60           1.047254   0.155899   6.718 1.85e-11 ***
## grupo_edad61+             1.053590   0.160399   6.569 5.08e-11 ***
## razaMex-Amer             -0.495218   0.143571  -3.449 0.000562 ***
## razaHispano              -0.249068   0.153928  -1.618 0.105645    
## razaNegro                -0.348555   0.108872  -3.202 0.001367 ** 
## razaAsiático             -0.655818   0.152987  -4.287 1.81e-05 ***
## razaOtra                 -0.106073   0.178135  -0.595 0.551533    
## educacionEdSecSin9       -0.279556   0.213060  -1.312 0.189486    
## educacionEdSecSinDip     -0.203128   0.166306  -1.221 0.221929    
## educacionEdSecDip        -0.108920   0.130534  -0.834 0.404045    
## educacionEdUnivInc        0.071669   0.116216   0.617 0.537444    
## estado_civilViudo/a      -0.051100   0.159164  -0.321 0.748169    
## estado_civilDivorciado/a  0.086586   0.128154   0.676 0.499269    
## estado_civilSeparado/a    0.386652   0.218560   1.769 0.076879 .  
## estado_civilSoltero/a     0.208309   0.131680   1.582 0.113663    
## estado_civilUnión libre   0.028427   0.160208   0.177 0.859164    
## quintil_pobrezaBajo       0.034519   0.140361   0.246 0.805738    
## quintil_pobrezaMedio     -0.141313   0.142583  -0.991 0.321641    
## quintil_pobrezaAlto      -0.129496   0.145299  -0.891 0.372802    
## quintil_pobrezaMuy alto   0.174312   0.149832   1.163 0.244672    
## actividad_vigorosaSí      0.004229   0.105669   0.040 0.968080    
## alcohol_12mSí             0.359024   0.173322   2.071 0.038319 *  
## fumo100Sí                 0.369749   0.085134   4.343 1.40e-05 ***
## imc                       0.036819   0.005682   6.480 9.16e-11 ***
## phq_sin_sueno             0.163138   0.011832  13.788  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4250.3  on 3509  degrees of freedom
## Residual deviance: 3744.8  on 3483  degrees of freedom
## AIC: 3798.8
## 
## Number of Fisher Scoring iterations: 4

library(sjPlot)

tab_model(modelo_log, transform = "exp", show.ci = 0.95,
  show.p = TRUE, show.se = TRUE, dv.labels = "Problema de sueño")

Para construir el modelo predictivo se utilizaron como predictores las variables cuantitativas edad, IMC, PAS1, HbA1c y puntuación total del PHQ-9, junto con las variables categóricas sexo, estado civil, quintil de pobreza, actividad vigorosa y antecedente de tabaquismo. La variable problema_sueno_consulta se definió como desenlace. La base de datos se dividió en una muestra de entrenamiento (70 %) para ajustar el modelo y seleccionar el punto de corte óptimo mediante el índice de Youden, y una muestra de validación independiente (30 %) para evaluar su desempeño predictivo fuera de la muestra, mediante AUC y matriz de confusión.

library(caret)
library(pROC)
library(sjPlot)

set.seed(1256)

indice_construccion <- createDataPartition(base_modelos$problema_sueno_consulta,p = 0.7,list = FALSE)

muestra_construccion <- base_modelos[indice_construccion, ]
muestra_validacion   <- base_modelos[-indice_construccion, ]

cat("Observaciones en muestra de construcción:", nrow(muestra_construccion), "\n")

## Observaciones en muestra de construcción: 2458

cat("Observaciones en muestra de validación  :", nrow(muestra_validacion), "\n")

## Observaciones en muestra de validación  : 1052

Revisamos la prevalencia de la variable de desenlace en cada base

cat("Prevalencia 'Sí' en muestra de construcción:",
    round(mean(muestra_construccion$problema_sueno_consulta == "Sí"), 4), "\n")

## Prevalencia 'Sí' en muestra de construcción: 0.2937

cat("Prevalencia 'Sí' en muestra de validación  :",
    round(mean(muestra_validacion$problema_sueno_consulta == "Sí"), 4), "\n")

## Prevalencia 'Sí' en muestra de validación  : 0.2937

modelo_completo <- glm(problema_sueno_consulta ~ grupo_edad + raza + educacion + estado_civil + 
                         quintil_pobreza + actividad_vigorosa + alcohol_12m + fumo100 + imc + phq_sin_sueno,
                       data = muestra_construccion, family = binomial)

modelo_step <- step(modelo_completo, direction = "both", trace = TRUE)

## Start:  AIC=2692.83
## problema_sueno_consulta ~ grupo_edad + raza + educacion + estado_civil + 
##     quintil_pobreza + actividad_vigorosa + alcohol_12m + fumo100 + 
##     imc + phq_sin_sueno
## 
##                      Df Deviance    AIC
## - estado_civil        5   2646.4 2690.4
## - quintil_pobreza     4   2644.5 2690.5
## - actividad_vigorosa  1   2638.9 2690.9
## - educacion           4   2645.8 2691.8
## <none>                    2638.8 2692.8
## - alcohol_12m         1   2642.9 2694.9
## - raza                5   2661.7 2705.7
## - fumo100             1   2656.8 2708.8
## - imc                 1   2670.5 2722.5
## - grupo_edad          3   2677.0 2725.0
## - phq_sin_sueno       1   2755.2 2807.2
## 
## Step:  AIC=2690.45
## problema_sueno_consulta ~ grupo_edad + raza + educacion + quintil_pobreza + 
##     actividad_vigorosa + alcohol_12m + fumo100 + imc + phq_sin_sueno
## 
##                      Df Deviance    AIC
## - quintil_pobreza     4   2652.1 2688.1
## - actividad_vigorosa  1   2646.6 2688.6
## - educacion           4   2653.9 2689.9
## <none>                    2646.4 2690.4
## - alcohol_12m         1   2650.2 2692.2
## + estado_civil        5   2638.8 2692.8
## - raza                5   2668.7 2702.7
## - fumo100             1   2664.8 2706.8
## - imc                 1   2677.3 2719.3
## - grupo_edad          3   2691.5 2729.5
## - phq_sin_sueno       1   2768.4 2810.4
## 
## Step:  AIC=2688.06
## problema_sueno_consulta ~ grupo_edad + raza + educacion + actividad_vigorosa + 
##     alcohol_12m + fumo100 + imc + phq_sin_sueno
## 
##                      Df Deviance    AIC
## - actividad_vigorosa  1   2652.2 2686.2
## <none>                    2652.1 2688.1
## - educacion           4   2660.2 2688.2
## - alcohol_12m         1   2655.9 2689.9
## + quintil_pobreza     4   2646.4 2690.4
## + estado_civil        5   2644.5 2690.5
## - raza                5   2674.8 2700.8
## - fumo100             1   2670.6 2704.6
## - imc                 1   2682.3 2716.3
## - grupo_edad          3   2698.8 2728.8
## - phq_sin_sueno       1   2777.2 2811.2
## 
## Step:  AIC=2686.17
## problema_sueno_consulta ~ grupo_edad + raza + educacion + alcohol_12m + 
##     fumo100 + imc + phq_sin_sueno
## 
##                      Df Deviance    AIC
## <none>                    2652.2 2686.2
## - educacion           4   2661.0 2687.0
## + actividad_vigorosa  1   2652.1 2688.1
## - alcohol_12m         1   2656.1 2688.1
## + quintil_pobreza     4   2646.6 2688.6
## + estado_civil        5   2644.6 2688.6
## - raza                5   2674.8 2698.8
## - fumo100             1   2670.6 2702.6
## - imc                 1   2682.3 2714.3
## - grupo_edad          3   2701.2 2729.2
## - phq_sin_sueno       1   2777.7 2809.7

summary(modelo_step)

## 
## Call:
## glm(formula = problema_sueno_consulta ~ grupo_edad + raza + educacion + 
##     alcohol_12m + fumo100 + imc + phq_sin_sueno, family = binomial, 
##     data = muestra_construccion)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -3.281724   0.331130  -9.911  < 2e-16 ***
## grupo_edad31-45       0.443277   0.171550   2.584 0.009768 ** 
## grupo_edad46-60       0.923697   0.164708   5.608 2.05e-08 ***
## grupo_edad61+         0.927044   0.160716   5.768 8.01e-09 ***
## razaMex-Amer         -0.501181   0.171653  -2.920 0.003503 ** 
## razaHispano          -0.070069   0.181156  -0.387 0.698912    
## razaNegro            -0.266299   0.126548  -2.104 0.035350 *  
## razaAsiático         -0.694296   0.184901  -3.755 0.000173 ***
## razaOtra              0.021715   0.205215   0.106 0.915729    
## educacionEdSecSin9   -0.503678   0.243708  -2.067 0.038759 *  
## educacionEdSecSinDip -0.341766   0.186078  -1.837 0.066256 .  
## educacionEdSecDip    -0.289905   0.143291  -2.023 0.043053 *  
## educacionEdUnivInc   -0.072243   0.131109  -0.551 0.581623    
## alcohol_12mSí         0.398727   0.206473   1.931 0.053467 .  
## fumo100Sí             0.432271   0.100773   4.290 1.79e-05 ***
## imc                   0.036253   0.006609   5.485 4.13e-08 ***
## phq_sin_sueno         0.140717   0.013094  10.747  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2976.5  on 2457  degrees of freedom
## Residual deviance: 2652.2  on 2441  degrees of freedom
## AIC: 2686.2
## 
## Number of Fisher Scoring iterations: 4

tab_model(modelo_completo, modelo_step, transform = NULL, show.ci = 0.95,
          show.p = TRUE, show.se = TRUE, 
          title = "Comparación de modelos de regresión logística (muestra de construcción)",
  dv.labels = c("Modelo completo", "Modelo Reducido (stepwise)"))

Primero en la muestra de construcción, luego en la muestra de validación

Predicción con umbral 0.5 en muestra de construcción

pred_log_construccion <- predict(modelo_step, newdata = muestra_construccion, type = "response")

muestra_construccion$pred_log <- ifelse(pred_log_construccion > 0.5, "Sí", "No")
muestra_construccion$pred_log <- factor(muestra_construccion$pred_log, levels = c("No", "Sí"))

tabla_logpred_construccion <- table(Predicho = muestra_construccion$pred_log,
  Real     = muestra_construccion$problema_sueno_consulta
)

confusionMatrix(tabla_logpred_construccion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho   No   Sí
##       No 1637  536
##       Sí   99  186
##                                           
##                Accuracy : 0.7417          
##                  95% CI : (0.7239, 0.7589)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 5.322e-05       
##                                           
##                   Kappa : 0.2437          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.25762         
##             Specificity : 0.94297         
##          Pos Pred Value : 0.65263         
##          Neg Pred Value : 0.75334         
##              Prevalence : 0.29373         
##          Detection Rate : 0.07567         
##    Detection Prevalence : 0.11595         
##       Balanced Accuracy : 0.60030         
##                                           
##        'Positive' Class : Sí              
## 

roc_obj_validacion1 <- roc(response  = muestra_validacion$problema_sueno_consulta,
                           predictor = pred_log_val,  levels = c("No", "Sí"))
auc_val <- auc(roc_obj_validacion1)
plot(roc_obj_validacion1,  col = "royalblue", lwd = 2, 
     main = "Curva ROC - Modelo Regresión Logística (muestra de validación)", 
     xlim = c(1.0, 0.0), ylim = c(0, 1), legacy.axes = TRUE, xaxs = "i", asp = NA)
legend("bottom", legend = paste0("AUC = ", round(auc_val, 3)),bty = "n")

Es una técnica estadística que se utiliza para clasificar individuos en grupos predefinidos basándose en varias variables predictoras. El objetivo es maximizar la separación entre grupos mientras se minimiza el error de clasificación.

Objetivos del Análisis Discriminante

1. Objetivo Descriptivo: Entender cómo las variables clasificadoras separan los grupos y cuál es la contribución de cada variable a esta separación.

2. Objetivo Predictivo: Establecer un criterio para asignar nuevos individuos a uno de los grupos basándose en sus características.

Para aplicar el Análisis Discriminante, los datos deben cumplir ciertos supuestos:

• Número de Individuos y Grupos: Se necesita un número adecuado de individuos en cada grupo. Idealmente, se recomienda tener al menos dos individuos por grupo.

• Variables: Las variables deben ser medidas en escala de intervalo o razón, no deben presentar multicolinealidad y deben seguir una distribución normal multivariante.

• Número de Variables: Para identificar el número adecuado de variables discriminantes, se debe cumplir la condición \(p < N - 2\), donde \(p\) es el número de variables y \(N\) es el número total de individuos.

Aunque la selección previa de variables puede realizarse comparando medias mediante pruebas t o ANOVA, se recomienda incluir todas las variables disponibles. De este modo, el Análisis Discriminante Lineal evalúa simultáneamente el poder predictor de todo el conjunto de datos sobre la variable de clasificación.

La regla discriminante se basa en la probabilidad a posteriori de pertenencia a un grupo. A cada individuo, se le calcula la probabilidad de pertenencia a cada grupo usando el Teorema de Bayes:

\[ P(G_m \mid x_i) = \frac{P(G_m) \, f(x_i \mid G_m)}{\sum_{k=1}^K P(G_k) \, f(x_i \mid G_k)} \]

donde:

• \(P(G_m \mid x_i)\): probabilidad posterior de que el individuo \(x_i\) pertenezca al grupo \(G_m\).

• \(P(G_m)\): probabilidad a priori del grupo \(G_m\).

• \(f(x_i \mid G_m)\): función de densidad de probabilidad de las variables predictoras \(x_i\), dado que el individuo pertenece al grupo \(G_m\).

• \(K\): número total de grupos o clases.

• El denominador \(\sum_{k=1}^K P(G_k) f(x_i \mid G_k)\) asegura la normalización para que las probabilidades posteriores sumen 1 entre todos los grupos.

• \(k\) es el número total de grupos.

El individuo se clasifica en el grupo con la mayor probabilidad a posteriori.

En el análisis discriminante, se pueden cometer dos tipos de errores:

1. Error Tipo I: Clasificar incorrectamente a un individuo en un grupo al que no pertenece.

2. Error Tipo II: No clasificar a un individuo en el grupo al que realmente pertenece.

Si se conocen los costes asociados a estos errores, se puede ajustar la regla discriminante ponderando las verosimilitudes por el coste de los errores.

Dos Grupos y una Variable Clasificadora

Para dos grupos y una variable \(X\), el objetivo es encontrar un punto de corte \(C\) que minimice el error de clasificación. La regla de clasificación se basa en la comparación del valor de \(X\) con \(C\). Si el valor de \(X\) es menor que \(C\), el individuo se clasifica en el grupo I, de lo contrario, en el grupo II.

Dos Grupos y Dos Variables Clasificadoras

Cuando se tienen dos variables \(X_1\) y \(X_2\), se busca una línea que separe los dos grupos en un espacio bidimensional. Se proyectan elipsoides de los grupos sobre un eje que combina ambas variables para minimizar el error de clasificación.

La función discriminante se puede expresar como:

\[ w_1 X_1 + w_2 X_2 - D = 0 \]

Donde \(w_1\) y \(w_2\) son los pesos de las variables \(X_1\) y \(X_2\), respectivamente, y \(D\) es el punto de corte.

Para \(p\) variables clasificadoras, la función discriminante se calcula como una combinación lineal de las variables:

\[ D = w_1 X_1 + w_2 X_2 + \cdots + w_p X_p \]

La puntuación discriminante para el individuo \(i\)-ésimo es:

\[ D_i = w_1 X_{1i} + w_2 X_{2i} + \cdots + w_p X_{pi} \]

El objetivo es maximizar la distancia entre los centros de los grupos y minimizar la varianza dentro de los grupos. La matriz de suma de cuadrados y productos cruzados (SCPC) se descompone en matrices entregrupos \(F\) y residual \(U\). La maximización se realiza usando:

\[ \frac{w' F w}{w' U w} \]

\(k\) Grupos y \(p\) Variables

En el caso de \(k\) grupos y \(p\) variables, se obtienen \(T = \min(k - 1, p)\) funciones discriminantes. La idea es encontrar combinaciones lineales de las variables que maximicen la separación entre los grupos. El criterio de maximización es similar al caso anterior:

\[ \frac{W' F W}{W' U W} \]

El LDA permite clasificar individuos según la presencia o ausencia de problemas de sueño.

A diferencia de la regresión logística, requiere variables cuantitativas con distribución aproximadamente normal multivariante y covarianzas similares entre los grupos. Se utilizan como predictores edad, IMC, PAS1, HbA1c y PHQ-total.

El modelo se ajusta con la muestra de construcción y se evalúa tanto en esta muestra como en la muestra de validación para valorar su desempeño fuera de la muestra.

library(MASS)
library(caret)

modelo_lda <- lda(problema_sueno_consulta ~ edad + imc + pas1 + hba1c + phq_sin_sueno,
  data = muestra_construccion)

modelo_lda

## Call:
## lda(problema_sueno_consulta ~ edad + imc + pas1 + hba1c + phq_sin_sueno, 
##     data = muestra_construccion)
## 
## Prior probabilities of groups:
##        No        Sí 
## 0.7062653 0.2937347 
## 
## Group means:
##        edad      imc     pas1    hba1c phq_sin_sueno
## No 49.17224 29.18007 125.4309 5.779781      2.040899
## Sí 53.81579 31.34529 128.5457 5.916205      4.139889
## 
## Coefficients of linear discriminants:
##                        LD1
## edad           0.026889637
## imc            0.059399222
## pas1          -0.001319246
## hba1c         -0.048202219
## phq_sin_sueno  0.233650761

Las probabilidades a priori reflejan la proporción de cada clase dentro de la muestra de construcción (deben coincidir con la prevalencia reportada más arriba, ~29.38% “Sí” y ~70.62% “No”, ya que el split fue estratificado):

modelo_lda$prior

##        No        Sí 
## 0.7062653 0.2937347

Perfil promedio por grupo (Group Means)

modelo_lda$means

##        edad      imc     pas1    hba1c phq_sin_sueno
## No 49.17224 29.18007 125.4309 5.779781      2.040899
## Sí 53.81579 31.34529 128.5457 5.916205      4.139889

Las medias de cada grupo describen las características biométricas y psicológicas promedio para la presencia (Sí) o ausencia (No) de problemas de sueño, calculadas dentro de la muestra de construcción.

modelo_lda$scaling

##                        LD1
## edad           0.026889637
## imc            0.059399222
## pas1          -0.001319246
## hba1c         -0.048202219
## phq_sin_sueno  0.233650761

Los coeficientes actúan como pesos aplicados a cada variable para construir la función discriminante lineal (LD1).

pred_lda_construccion <- predict(modelo_lda)  # sin newdata = usa muestra_construccion
tabla_lda_construccion <- table(
  Predicho = pred_lda_construccion$class,
  Real     = muestra_construccion$problema_sueno_consulta)
confusionMatrix(tabla_lda_construccion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho   No   Sí
##       No 1629  565
##       Sí  107  157
##                                           
##                Accuracy : 0.7266          
##                  95% CI : (0.7085, 0.7442)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 0.01376         
##                                           
##                   Kappa : 0.1912          
##                                           
##  Mcnemar's Test P-Value : < 2e-16         
##                                           
##             Sensitivity : 0.21745         
##             Specificity : 0.93836         
##          Pos Pred Value : 0.59470         
##          Neg Pred Value : 0.74248         
##              Prevalence : 0.29373         
##          Detection Rate : 0.06387         
##    Detection Prevalence : 0.10740         
##       Balanced Accuracy : 0.57791         
##                                           
##        'Positive' Class : Sí              
## 

pred_lda_validacion <- predict(modelo_lda, newdata = muestra_validacion)
tabla_lda_validacion <- table(
  Predicho = pred_lda_validacion$class,
  Real     = muestra_validacion$problema_sueno_consulta)
confusionMatrix(tabla_lda_validacion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho  No  Sí
##       No 714 248
##       Sí  29  61
##                                          
##                Accuracy : 0.7367         
##                  95% CI : (0.709, 0.7631)
##     No Information Rate : 0.7063         
##     P-Value [Acc > NIR] : 0.01578        
##                                          
##                   Kappa : 0.1997         
##                                          
##  Mcnemar's Test P-Value : < 2e-16        
##                                          
##             Sensitivity : 0.19741        
##             Specificity : 0.96097        
##          Pos Pred Value : 0.67778        
##          Neg Pred Value : 0.74220        
##              Prevalence : 0.29373        
##          Detection Rate : 0.05798        
##    Detection Prevalence : 0.08555        
##       Balanced Accuracy : 0.57919        
##                                          
##        'Positive' Class : Sí             
## 

Por defecto, predict.lda() clasifica usando un corte implícito de 0.5 sobre la probabilidad posterior. Igual que en la regresión logística, buscamos el punto de corte óptimo con la curva ROC de la muestra de construcción y lo aplicamos, a la muestra de validación.

posterior_si_construccion <- pred_lda_construccion$posterior[, "Sí"]
roc_obj_lda_construccion <- roc(
  response  = muestra_construccion$problema_sueno_consulta,
  predictor = posterior_si_construccion,
  levels    = c("No", "Sí"))

plot(roc_obj_lda_construccion, col  = "seagreen", lwd  = 2,
     main = "Curva ROC - LDA (muestra de construcción)")

corte_optimo_lda <- coords(roc_obj_lda_construccion, x = "best",
                           best.method = "youden", 
                           ret = c("threshold", "sensitivity", "specificity"))
print(corte_optimo_lda)

##   threshold sensitivity specificity
## 1  0.273847   0.6177285   0.6906682

umbral_lda <- corte_optimo_lda$threshold

Clasificación en la muestra de construcción con el umbral óptimo

muestra_construccion$pred_lda_opt <- ifelse(posterior_si_construccion > umbral_lda, "Sí", "No")
muestra_construccion$pred_lda_opt <- factor(muestra_construccion$pred_lda_opt, levels = levels(muestra_construccion$problema_sueno_consulta))

tabla_lda_opt_construccion <- table(
  Predicho = muestra_construccion$pred_lda_opt,
  Real     = muestra_construccion$problema_sueno_consulta)

confusionMatrix(tabla_lda_opt_construccion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho   No   Sí
##       No 1199  276
##       Sí  537  446
##                                           
##                Accuracy : 0.6692          
##                  95% CI : (0.6502, 0.6878)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.2789          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.6177          
##             Specificity : 0.6907          
##          Pos Pred Value : 0.4537          
##          Neg Pred Value : 0.8129          
##              Prevalence : 0.2937          
##          Detection Rate : 0.1814          
##    Detection Prevalence : 0.3999          
##       Balanced Accuracy : 0.6542          
##                                           
##        'Positive' Class : Sí              
## 

posterior_si_validacion <- pred_lda_validacion$posterior[, "Sí"]

roc_obj_lda_validacion <- roc(
  response  = muestra_validacion$problema_sueno_consulta,
  predictor = posterior_si_validacion,
  levels    = c("No", "Sí"))

plot(roc_obj_lda_validacion, col  = "darkorange",
  lwd  = 2, main = "Curva ROC - LDA (muestra de validación)")

muestra_validacion$pred_lda_opt <- ifelse(posterior_si_validacion > umbral_lda, "Sí", "No")
muestra_validacion$pred_lda_opt <- factor(muestra_validacion$pred_lda_opt, levels = levels(muestra_construccion$problema_sueno_consulta))

tabla_lda_opt_validacion <- table(
  Predicho = muestra_validacion$pred_lda_opt,
  Real     = muestra_validacion$problema_sueno_consulta)

confusionMatrix(tabla_lda_opt_validacion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho  No  Sí
##       No 521 109
##       Sí 222 200
##                                           
##                Accuracy : 0.6854          
##                  95% CI : (0.6563, 0.7133)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 0.9354          
##                                           
##                   Kappa : 0.3148          
##                                           
##  Mcnemar's Test P-Value : 7.457e-10       
##                                           
##             Sensitivity : 0.6472          
##             Specificity : 0.7012          
##          Pos Pred Value : 0.4739          
##          Neg Pred Value : 0.8270          
##              Prevalence : 0.2937          
##          Detection Rate : 0.1901          
##    Detection Prevalence : 0.4011          
##       Balanced Accuracy : 0.6742          
##                                           
##        'Positive' Class : Sí              
## 

cat("AUC LDA construcción:", round(auc(roc_obj_lda_construccion), 4), "\n")

## AUC LDA construcción: 0.6956

cat("AUC LDA validación  :", round(auc(roc_obj_lda_validacion), 4), "\n")

## AUC LDA validación  : 0.7323

Aunque el modelo presenta una precisión general aceptable (73.94%), su rendimiento real esta bastante desequilibrado, ya que muestra un alto poder para descartar correctamente a los casos negativos con una especificidad del 92.99% y un valor predictivo negativo del 75.68%, pero una deficiente capacidad para identificar a los individuos con la condición, donde apenas detecta un 28.15% de los casos verdaderos positivos (sensibilidad).

Además, se observa un bajo por no decir débil índice de concordancia ajustado por azar (Kappa = 0.2519) y una exactitud balanceada de solo 60.57%, confirman que el modelo tiende a sobrepredecir la ausencia de la condición.

Los árboles de clasificación y regresión son métodos no paramétricos utilizados para clasificar, discriminar o predecir unidades de observación a partir de un conjunto de variables predictoras. Son técnicas de dependencia, en las que se busca explicar o predecir una variable respuesta (Y) a partir de un conjunto de variables independientes (X_1,,X_p).

La elección del tipo de árbol depende de la naturaleza de la variable respuesta:

• Árboles de regresión: se utilizan cuando \(Y\) es cuantitativa y buscan estimar el valor esperado de la respuesta:

\[ E(Y\mid X_1,\ldots,X_p) \]

• Árboles de clasificación: se utilizan cuando \(Y\) es categórica y buscan estimar la probabilidad de pertenecer a cada categoría:

\[ P(Y=y\mid X_1,\ldots,X_p) \]

Entre sus principales ventajas, requieren pocos supuestos distribucionales, pueden capturar relaciones no lineales e interacciones entre variables y son relativamente robustos frente a valores atípicos.

El procedimiento se conoce como partición binaria recursiva. Inicialmente, todas las observaciones se encuentran en un único nodo o grupo. El algoritmo identifica la división que permite separar mejor las observaciones según la variable respuesta y divide el nodo en dos grupos. Este proceso se repite recursivamente en los nuevos nodos hasta alcanzar un criterio de parada, generando una estructura en forma de árbol que permite realizar predicciones y caracterizar los grupos obtenidos.

Un árbol de clasificación o regresión está compuesto por diferentes tipos de nodos:

• Nodo raíz: corresponde al nodo inicial que contiene todas las observaciones.
• Nodos internos: representan puntos donde se realiza una nueva partición de las observaciones.
• Nodos terminales: son los nodos finales, que ya no se dividen. Se busca que sean lo más homogéneos posible respecto a la variable respuesta, utilizando medidas de impureza.

El algoritmo busca dividir cada nodo en dos grupos mutuamente excluyentes, utilizando una variable predictora y un punto de corte. Por ejemplo, una variable como edad podría dividirse en personas de 45 años o más y personas menores de 45 años.

En cada partición se selecciona una única variable y el punto de corte que produzcan la mayor reducción de la impureza, generando grupos lo más homogéneos posible.

El procedimiento es recursivo porque cada nodo resultante puede volver a dividirse mediante nuevas variables y puntos de corte. El proceso continúa hasta alcanzar algún criterio de parada, como un tamaño mínimo de nodo o un nivel determinado de complejidad.

Existen diferentes medidas para evaluar la calidad de las particiones. En árboles de clasificación, una de las más utilizadas es el índice de Gini, que cuantifica la impureza de las categorías dentro de cada nodo.

Mide la impureza de un nodo, es decir, qué tan mezcladas se encuentran las diferentes clases dentro de una partición. Supongamos que tenemos dos clases: Clase 1 y Clase 2.

\[ IG=1-p_1^2-p_2^2 \]

Donde \(p_1\) es la proporción de observaciones que pertenecen a la Clase 1 y \(p_2\) es la proporción de observaciones que pertenecen a la Clase 2.

El índice de Gini alcanza su mayor valor cuando las clases están igualmente representadas. Para dos clases:

\[ IG=1-0.5^2-0.5^2=0.5 \]

Por el contrario, cuando todas las observaciones pertenecen a una sola clase, el nodo es completamente puro y el índice de Gini es igual a cero:

\[ IG=1-1^2-0^2=0 \]

Por ejemplo, si en una partición inicial hay 296 mujeres y 61 hombres, para un total de 357 observaciones, el índice de Gini es:

\[ IG=1-(296/357)^2-(61/357)^2=0.2834 \]

Posteriormente, se consideran todos los posibles cortes para todas las variables predictoras y se selecciona aquel que produzca la mayor reducción de la impureza, es decir, el que genere los nodos resultantes más homogéneos.

El índice de Gini después de una división se obtiene como un promedio ponderado de la impureza de los nodos resultantes, donde las ponderaciones corresponden a la proporción de observaciones que queda en cada nodo.

Por ejemplo, si una variable como sexo genera directamente dos grupos, se calcula el Gini de cada grupo y posteriormente se obtiene el Gini ponderado de la división. Si la variable es cuantitativa, como años de escolaridad, se evalúan diferentes puntos de corte, por ejemplo, más de 10 años frente a 10 años o menos.

Para \(K\) categorías, el índice de Gini se expresa como:

\[ IG = 1-\sum_{k=1}^K p_k^2 \]

• Valores cercanos a 0: indican nodos más puros y homogéneos.
• Valores mayores: indican mayor mezcla entre las categorías.
• Una clasificación perfectamente homogénea presenta un índice de Gini igual a 0.

Es otra medida de impureza utilizada para evaluar la calidad de las particiones. Mide el grado de incertidumbre existente respecto a la categoría de la variable respuesta.

\[ H=-\sum_{k=1}^K p_k\log_2(p_k) \]

Una entropía baja indica que las observaciones están concentradas principalmente en una categoría, mientras que una entropía alta indica mayor mezcla entre las categorías.

La ganancia de información corresponde a la reducción de la entropía producida por una división:

\[ \text{Ganancia de información} = \text{Entropía del nodo padre} - \text{Entropía ponderada de los nodos hijos} \]

Por tanto, se prefieren las particiones que produzcan una mayor ganancia de información, ya que generan nodos con menor incertidumbre.

En este análisis, las variables que originalmente son continuas se mantendrán como variables continuas, permitiendo que el algoritmo determine automáticamente los puntos de corte más adecuados.

Al igual que en los modelos anteriores, el árbol se ajusta únicamente utilizando la muestra de construcción (70%). Posteriormente, su desempeño se evalúa tanto en la muestra de construcción como, de manera independiente, en la muestra de validación (30%).

Esta estrategia permite evaluar el desempeño en la muestra de construcción y, principalmente, determinar qué tan bien el árbol generaliza a datos que no fueron utilizados durante su entrenamiento.

modelo_cart <- rpart(problema_sueno_consulta ~ edad + sexo + raza + estado_civil + 
                       quintil_pobreza + imc + pas1 + hba1c + actividad_vigorosa + 
                       alcohol_12m + fumo100 + phq_sin_sueno, 
                     data = muestra_construccion, method = "class")

rpart.plot(modelo_cart, type = 2, fallen.leaves = TRUE, cex = 0.7,
           main = "Árbol CART sin podar (muestra de construcción)")

modelo_cart$variable.importance

## phq_sin_sueno          edad           imc         hba1c          pas1 
##   62.94466907    9.43313805    7.30535060    0.88238309    0.15214739 
##   alcohol_12m 
##    0.03844921

printcp(modelo_cart)

## 
## Classification tree:
## rpart(formula = problema_sueno_consulta ~ edad + sexo + raza + 
##     estado_civil + quintil_pobreza + imc + pas1 + hba1c + actividad_vigorosa + 
##     alcohol_12m + fumo100 + phq_sin_sueno, data = muestra_construccion, 
##     method = "class")
## 
## Variables actually used in tree construction:
## [1] edad          imc           phq_sin_sueno
## 
## Root node error: 722/2458 = 0.29373
## 
## n= 2458 
## 
##         CP nsplit rel error  xerror     xstd
## 1 0.022853      0   1.00000 1.00000 0.031276
## 2 0.010000      4   0.89889 0.95014 0.030801

plotcp(modelo_cart)

se elige el cp que minimiza el error de validación cruzada interna (xerror) y se poda el árbol con prune().

cptable <- modelo_cart$cptable

cp_optimo <- cptable[which.min(cptable[, "xerror"]), "CP"]

# min_xerror <- min(cptable[, "xerror"])
# se_min     <- cptable[which.min(cptable[, "xerror"]), "xstd"]
# cp_optimo  <- cptable[cptable[, "xerror"] <= min_xerror + se_min, "CP"][1]

cat("cp óptimo seleccionado:", cp_optimo, "\n")

## cp óptimo seleccionado: 0.01

modelo_cart_podado <- prune(modelo_cart, cp = cp_optimo)

rpart.plot(modelo_cart_podado, type = 2, fallen.leaves = TRUE, cex = 0.7,
           main = "Árbol CART podado (muestra de construcción)")

modelo_cart_podado$variable.importance

## phq_sin_sueno          edad           imc         hba1c          pas1 
##   62.94466907    9.43313805    7.30535060    0.88238309    0.15214739 
##   alcohol_12m 
##    0.03844921

pred_cart_construccion <- predict(modelo_cart_podado, newdata = muestra_construccion, type = "class")
tabla_cart_construccion <- table(Predicho = pred_cart_construccion,
  Real     = muestra_construccion$problema_sueno_consulta)
confusionMatrix(tabla_cart_construccion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho   No   Sí
##       No 1628  541
##       Sí  108  181
##                                           
##                Accuracy : 0.736           
##                  95% CI : (0.7181, 0.7533)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 0.000596        
##                                           
##                   Kappa : 0.2285          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.25069         
##             Specificity : 0.93779         
##          Pos Pred Value : 0.62630         
##          Neg Pred Value : 0.75058         
##              Prevalence : 0.29373         
##          Detection Rate : 0.07364         
##    Detection Prevalence : 0.11758         
##       Balanced Accuracy : 0.59424         
##                                           
##        'Positive' Class : Sí              
## 

pred_cart_validacion <- predict(modelo_cart_podado, newdata = muestra_validacion, type = "class")
tabla_cart_validacion <- table(Predicho = pred_cart_validacion, 
                               Real     = muestra_validacion$problema_sueno_consulta)
confusionMatrix(tabla_cart_validacion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho  No  Sí
##       No 710 238
##       Sí  33  71
##                                           
##                Accuracy : 0.7424          
##                  95% CI : (0.7148, 0.7686)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 0.005147        
##                                           
##                   Kappa : 0.2299          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.22977         
##             Specificity : 0.95559         
##          Pos Pred Value : 0.68269         
##          Neg Pred Value : 0.74895         
##              Prevalence : 0.29373         
##          Detection Rate : 0.06749         
##    Detection Prevalence : 0.09886         
##       Balanced Accuracy : 0.59268         
##                                           
##        'Positive' Class : Sí              
## 

Igual que en regresión logística y en LDA, type = "class" clasifica con un corte implícito de 0.5 sobre la probabilidad.

prob_cart_construccion <- predict(modelo_cart_podado, newdata=muestra_construccion, type = "prob")[,"Sí"]
prob_cart_validacion   <- predict(modelo_cart_podado, newdata=muestra_validacion,   type = "prob")[,"Sí"]

roc_cart_construccion <- roc(
  response  = muestra_construccion$problema_sueno_consulta,
  predictor = prob_cart_construccion,
  levels    = c("No", "Sí"))

plot(roc_cart_construccion, col = "seagreen", lwd = 2,
     main = "Curva ROC - CART (muestra de construcción)")

corte_optimo_cart <- coords(roc_cart_construccion, x = "best", best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"))

print(corte_optimo_cart)

##   threshold sensitivity specificity
## 1 0.3134965   0.3434903   0.8761521

umbral_cart <- corte_optimo_cart$threshold

roc_cart_validacion <- roc(response  = muestra_validacion$problema_sueno_consulta,
  predictor = prob_cart_validacion,levels    = c("No", "Sí"))

plot(roc_cart_construccion, col = "seagreen", lwd = 2,
     main = "Curva ROC - CART (construcción vs. validación)")
plot(roc_cart_validacion, col = "darkorange", lwd = 2, add = TRUE)
legend("bottomright", legend = c("Construcción", "Validación"),
       col = c("seagreen", "darkorange"), lwd = 2)

cat("AUC CART construcción:", round(auc(roc_cart_construccion), 4), "\n")

## AUC CART construcción: 0.6152

cat("AUC CART validación  :", round(auc(roc_cart_validacion), 4), "\n")

## AUC CART validación  : 0.6299

Son modelos computacionales utilizados para aprender relaciones complejas entre un conjunto de variables de entrada y una variable de respuesta. Se inspiran en el funcionamiento de las redes neuronales biológicas.

Son métodos de aprendizaje supervisado que pueden utilizarse tanto para clasificación como para regresión. Su principal característica es que pueden representar relaciones no lineales y complejas entre las variables predictoras y la respuesta.

Una red neuronal puede entenderse como una sucesión de operaciones matemáticas en las que las variables de entrada se combinan mediante pesos y se transforman a través de funciones de activación. Durante el entrenamiento, estos pesos se ajustan para minimizar el error de las predicciones.

• Modelar relaciones complejas y no lineales entre variables.
• Realizar predicciones para problemas de clasificación o regresión.
• Identificar patrones difíciles de representar mediante modelos estadísticos tradicionales.
• Combinar múltiples variables predictoras en un único modelo.
• Obtener predicciones precisas en problemas con estructuras complejas.
• Servir como herramienta de aprendizaje automático para problemas de alta dimensionalidad.

Una red neuronal está formada por diferentes capas de unidades o “neuronas”:

• Capa de entrada: recibe las variables predictoras.
• Capas ocultas: transforman progresivamente la información mediante combinaciones de las variables de entrada.
• Capa de salida: produce la predicción final del modelo.

De forma esquemática:

\[ \text{Entradas} \rightarrow \text{Capas ocultas} \rightarrow \text{Salida} \]

Cada conexión entre neuronas tiene asociado un peso, que determina la influencia de una entrada sobre la siguiente unidad. Una neurona recibe varias entradas, las combina mediante una suma ponderada y posteriormente aplica una función de activación:

\[ z = \sum_{j=1}^{p} w_jx_j+b \]

\[ a = f(z) \]

donde \(x_j\) representa una variable de entrada, \(w_j\) su peso, \(b\) el sesgo (bias) y \(f(\cdot)\) la función de activación.

Las funciones de activación introducen no linealidad en la red, permitiendo que esta represente relaciones más complejas que una combinación lineal de las variables. Entre las funciones más utilizadas se encuentran:

Sigmoide

\[ f(z)=\frac{1}{1+e^{-z}} \]

Produce valores entre 0 y 1 y es especialmente utilizada en problemas de clasificación binaria.

ReLU (Unidad Lineal Rectificada)

\[ f(z)=\max(0,z) \]

Es una función de activación para capas ocultas que convierte los valores negativos en cero y mantiene los positivos sin cambios.

Softmax

Función de activación usada en la capa de salida para problemas de clasificación multiclase. Transforma los valores de la red en una distribución de probabilidades que suman 1 (o 100%)

Proceso para ajustar pesos y sesgos con el fin de minimizar el error. Evalúa la precisión comparando las predicciones con los valores reales mediante una función de pérdida (loss function).

Por ejemplo, en clasificación binaria puede utilizarse la entropía cruzada:

\[ L = -\left[ y\log(\hat{p})+ (1-y)\log(1-\hat{p}) \right] \]

donde \(y\) es el valor observado y \(\hat{p}\) la probabilidad predicha por la red.

El proceso de entrenamiento busca minimizar esta función de pérdida.

Durante la propagación hacia adelante (forward propagation), los datos pasan desde la capa de entrada hasta la capa de salida, generando una predicción. Posteriormente se calcula el error entre la predicción y el valor observado.

Mediante la retropropagación del error (backpropagation), se calcula cómo debe modificarse cada peso para reducir el error. Los pesos se actualizan mediante algoritmos de optimización, como el descenso del gradiente:

\[ w^{(t+1)} = w^{(t)} - \eta \frac{\partial L}{\partial w} \]

donde \(\eta\) es la tasa de aprendizaje y \(\frac{\partial L}{\partial w}\) representa el gradiente de la función de pérdida respecto al peso.

Este proceso se repite durante múltiples iteraciones hasta que el modelo alcanza un nivel de error adecuado.

La capacidad de una red neuronal depende, entre otros aspectos, de su arquitectura.

Algunos elementos importantes son:

• Número de capas ocultas.
• Número de neuronas por capa.
• Funciones de activación.
• Tasa de aprendizaje.
• Número de iteraciones o épocas (epochs).
• Tamaño de los lotes (batch size).
• Función de pérdida.

Una red demasiado sencilla puede ser incapaz de representar adecuadamente la relación entre las variables. Por el contrario, una red excesivamente compleja puede aprender demasiado bien los datos de entrenamiento y presentar un rendimiento deficiente sobre nuevos datos.

Uno de los principales problemas de las redes neuronales es el sobreajuste (overfitting), el cual ocurre cuando la red aprende características específicas del conjunto de entrenamiento, incluyendo ruido o patrones particulares, pero pierde capacidad para generalizar a nuevos individuos.

Para mitigar este problema y lograr una alta capacidad de generalización, se aplican técnicas como la división en conjuntos (entrenamiento, validación y prueba), la validación cruzada, la simplificación de la arquitectura, la detención temprana, la regularización de pesos y el dropout.

Las redes neuronales también son sensibles a la escala de las variables numéricas. Por ello, es habitual estandarizar o normalizar los predictores antes del entrenamiento. Una transformación frecuente es la estandarización:

\[ z_i=\frac{x_i-\bar{x}}{s_x} \]

donde \(\bar{x}\) es la media y \(s_x\) la desviación estándar de la variable.

Las variables categóricas también deben transformarse a una representación numérica adecuada, por ejemplo mediante variables indicadoras. Es importante realizar estas transformaciones utilizando únicamente la información disponible en el conjunto de entrenamiento para evitar fuga de información.

El desempeño de una red neuronal se evalúa sobre datos no utilizados durante el ajuste de parámetros, eligiendo las métricas según el objetivo del análisis y las características del problema. En clasificación, suelen emplearse la exactitud, sensibilidad, especificidad, precisión, F1-score, el área bajo la curva ROC (AUC) y la entropía cruzada; mientras que en problemas de regresión se recurre a métricas como el error cuadrático medio (MSE), su raíz (RMSE), el error absoluto medio (MAE) y el coeficiente de determinación (\(R^2\)).

Al considerarse modelos de «caja negra», las redes neuronales no admiten una interpretación directa de sus parámetros como una regresión, ya que la información se distribuye entre múltiples capas, pesos y funciones de activación; por ello, la interpretación del modelo exige recurrir a:

• Permutation Importance: Mezcla los datos de una variable al azar; si la precisión se cae, la variable es relevante; si no cambia, es irrelevante.

• Valores SHAP: Miden cuánto suma o resta cada variable en la predicción de un caso individual, evaluando su aporte exacto al resultado final.

• Gráficos de dependencia parcial: Muestran cómo varía la predicción a medida que sube o baja una variable específica, dejando las demás constantes.

• Explicación local: Aproxima el comportamiento de la red mediante un modelo simple (como una regresión) alrededor de una sola predicción para explicar ese caso concreto.

Las redes neuronales son especialmente útiles cuando el objetivo principal es la predicción y existe una estructura compleja o no lineal que otros modelos más simples no consiguen representar adecuadamente.

Usaremos exactamente los mismos predictores que en la regresión logística, para que la comparación entre métodos sea justa. Las redes neuronales necesitan variables numéricas estandarizadas y deben calcularse solo con la muestra de construcción y aplicarse tal cual a la de validación.

variables_nn <- c("edad", "sexo", "quintil_pobreza", "imc", "pas1",
  "hba1c", "actividad_vigorosa", "fumo100", "phq_total")

datos_nn_construccion <- muestra_construccion %>%
  mutate( sexo = as.numeric(sexo), quintil_pobreza = as.numeric(quintil_pobreza),
    actividad_vigorosa = as.numeric(actividad_vigorosa),     fumo100 = as.numeric(fumo100) )

datos_nn_validacion <- muestra_validacion %>%
  mutate(sexo = as.numeric(sexo), quintil_pobreza = as.numeric(quintil_pobreza),
    actividad_vigorosa = as.numeric(actividad_vigorosa), fumo100 = as.numeric(fumo100)  )

escalado_construccion <- scale(datos_nn_construccion[variables_nn])

centro <- attr(escalado_construccion, "scaled:center")
escala <- attr(escalado_construccion, "scaled:scale")

datos_nn_construccion[variables_nn] <- escalado_construccion
datos_nn_validacion[variables_nn]   <- scale(datos_nn_validacion[variables_nn],
                                             center = centro, scale  = escala)

set.seed(57123)
modelo_nn <- nnet(problema_sueno_consulta ~ edad + sexo + quintil_pobreza + imc + pas1 +
    hba1c + actividad_vigorosa + fumo100 + phq_total, data   = datos_nn_construccion,
  size   = 5, decay  = 0.01, maxit  = 500, trace  = FALSE)

modelo_nn

## a 9-5-1 network with 56 weights
## inputs: edad sexo quintil_pobreza imc pas1 hba1c actividad_vigorosa fumo100 phq_total 
## output(s): problema_sueno_consulta 
## options were - entropy fitting  decay=0.01

plotnet(modelo_nn, alpha = 0.6, cex.val = 0.8, cex.lab = 0.8,
        main = "Red neuronal artificial (muestra de construcción)")

pred_nn_construccion <- predict(modelo_nn, newdata = datos_nn_construccion, type = "class")
pred_nn_construccion <- factor(pred_nn_construccion, levels = levels(muestra_construccion$problema_sueno_consulta))
tabla_nn_construccion <- table(Predicho = pred_nn_construccion,
  Real     = datos_nn_construccion$problema_sueno_consulta)
confusionMatrix(tabla_nn_construccion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho   No   Sí
##       No 1566  445
##       Sí  170  277
##                                           
##                Accuracy : 0.7498          
##                  95% CI : (0.7322, 0.7668)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 8.457e-07       
##                                           
##                   Kappa : 0.3215          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.3837          
##             Specificity : 0.9021          
##          Pos Pred Value : 0.6197          
##          Neg Pred Value : 0.7787          
##              Prevalence : 0.2937          
##          Detection Rate : 0.1127          
##    Detection Prevalence : 0.1819          
##       Balanced Accuracy : 0.6429          
##                                           
##        'Positive' Class : Sí              
## 

pred_nn_validacion <- predict(modelo_nn, newdata = datos_nn_validacion, type = "class")
pred_nn_validacion <- factor(pred_nn_validacion, levels = levels(muestra_construccion$problema_sueno_consulta))
tabla_nn_validacion <- table(Predicho = pred_nn_validacion,
  Real     = datos_nn_validacion$problema_sueno_consulta)
confusionMatrix(tabla_nn_validacion, positive = "Sí")

## Confusion Matrix and Statistics
## 
##         Real
## Predicho  No  Sí
##       No 683 201
##       Sí  60 108
##                                           
##                Accuracy : 0.7519          
##                  95% CI : (0.7246, 0.7777)
##     No Information Rate : 0.7063          
##     P-Value [Acc > NIR] : 0.0005507       
##                                           
##                   Kappa : 0.3101          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.3495          
##             Specificity : 0.9192          
##          Pos Pred Value : 0.6429          
##          Neg Pred Value : 0.7726          
##              Prevalence : 0.2937          
##          Detection Rate : 0.1027          
##    Detection Prevalence : 0.1597          
##       Balanced Accuracy : 0.6344          
##                                           
##        'Positive' Class : Sí              
## 

En la muestra de construcción y punto de corte óptimo (Youden)

prob_nn_construccion <- predict(modelo_nn, newdata = datos_nn_construccion, type = "raw")[, 1]
prob_nn_validacion   <- predict(modelo_nn, newdata = datos_nn_validacion,   type = "raw")[, 1]

roc_nn_construccion <- roc(response  = datos_nn_construccion$problema_sueno_consulta,
  predictor = prob_nn_construccion, levels    = c("No", "Sí"))

plot(roc_nn_construccion, col = "seagreen", lwd = 2,
     main = "Curva ROC - Red neuronal (muestra de construcción)")

corte_optimo_nn <- coords(roc_nn_construccion, x = "best",
  best.method = "youden", ret = c("threshold", "sensitivity", "specificity"))

print(corte_optimo_nn)

##   threshold sensitivity specificity
## 1 0.2812331   0.7257618    0.687212

umbral_nn <- corte_optimo_nn$threshold

roc_nn_validacion <- roc(response  = datos_nn_validacion$problema_sueno_consulta,
  predictor = prob_nn_validacion, levels    = c("No", "Sí"))

plot(roc_nn_construccion, col = "seagreen", lwd = 2,
     main = "Curva ROC - Red neuronal (construcción vs. validación)")
plot(roc_nn_validacion, col = "darkorange", lwd = 2, add = TRUE)
legend("bottomright", legend = c("Construcción", "Validación"),
       col = c("seagreen", "darkorange"), lwd = 2)

cat("AUC red neuronal construcción:", round(auc(roc_nn_construccion), 4), "\n")

## AUC red neuronal construcción: 0.7696

cat("AUC red neuronal validación  :", round(auc(roc_nn_validacion), 4), "\n")

## AUC red neuronal validación  : 0.7507