Tema 3. Modelo Analítico

3. Modelo Analítico

library(readxl)       # leer el archivo de Excel
library(dplyr)        # manipular datos
library(tidyr)        # limpiar datos (drop_na)
library(ggplot2)      # gráficos
library(plotly)       # gráficos interactivos y 3D
library(skimr)        # resumen descriptivo completo
library(rsample)      # dividir en entrenamiento y prueba
library(broom)        # extraer coeficientes y residuos en tablas
library(yardstick)    # métricas de desempeño (AUC, sensibilidad, RMSE)
library(randomForest) # bosque aleatorio (algoritmo alternativo)
library(DiagrammeR)   # diagrama del esquema conceptual

sueno <- read_excel("Sleep_health_and_lifestyle_dataset.xlsx")
set.seed(2026)

3.1 Contexto y objetivos

3.1.1 Problema de gestión

Una IPS que presta servicios de salud ocupacional y promoción y prevención a empresas quiere identificar a los trabajadores con trastornos del sueño (insomnio o apnea). Estos trastornos aumentan el ausentismo, los accidentes laborales y el riesgo cardiovascular, pero confirmarlos requiere estudios costosos (polisomnografía) y valoración por especialista, con cupos limitados.

Reto: construir un modelo que, con datos que se recogen en un examen ocupacional de rutina (edad, sexo, presión arterial, IMC, nivel de estrés, horas de sueño, actividad física), estime la probabilidad de que un trabajador tenga un trastorno del sueño. Con esa probabilidad, la IPS puede priorizar a quién enviar primero a estudio especializado y dirigir mejor los recursos.

• Modelo principal (clasificación): regresión logística (técnica que estima la probabilidad de un resultado Sí/No) para predecir Sleep Disorder.
• Modelo complementario (regresión): regresión lineal para explicar la calidad del sueño (escala 1 a 10) e identificar qué factores modificables (estrés, actividad física) conviene intervenir en los programas de bienestar.

3.1.2 Objetivos y métricas de éxito

Modelo	Métrica	Qué mide	Meta
Clasificación	AUC (área bajo la curva ROC)	Capacidad de separar a quienes tienen y no tienen trastorno (0,5 = azar, 1 = perfecto)	≥ 0,85 en prueba
Clasificación	Sensibilidad	De los que sí tienen trastorno, cuántos detecta el modelo	≥ 0,80 (no dejar casos sin remitir)
Clasificación	Especificidad	De los que no tienen trastorno, cuántos descarta correctamente	≥ 0,75 (no desperdiciar cupos de estudio)
Regresión	RMSE (raíz del error cuadrático medio)	Error típico al predecir la calidad del sueño, en puntos	≤ 0,5 puntos
Ambos	Brecha entrenamiento–prueba	Si el modelo “memorizó” los datos (sobreajuste)	Diferencia pequeña (AUC ≤ 0,05)

En tamizaje, la sensibilidad pesa más que la especificidad: es preferible remitir a alguien sano a un estudio que dejar sin diagnóstico a un trabajador con apnea.

3.2 Selección de variables

3.2.1 Exploración de la base con glimpse() y skim()

glimpse(sueno)

## Rows: 374
## Columns: 13
## $ `Person ID`               <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1…
## $ Gender                    <chr> "Male", "Male", "Male", "Male", "Male", "Mal…
## $ Age                       <dbl> 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, …
## $ Occupation                <chr> "Software Engineer", "Doctor", "Doctor", "Sa…
## $ `Sleep Duration`          <chr> "6.1", "6.2", "6.2", "5.9", "5.9", "5.9", "6…
## $ `Quality of Sleep`        <dbl> 6, 6, 6, 4, 4, 4, 6, 7, 7, 7, 6, 7, 6, 6, 6,…
## $ `Physical Activity Level` <dbl> 42, 60, 60, 30, 30, 30, 40, 75, 75, 75, 30, …
## $ `Stress Level`            <dbl> 6, 8, 8, 8, 8, 8, 7, 6, 6, 6, 8, 6, 8, 8, 8,…
## $ `BMI Category`            <chr> "Overweight", "Normal", "Normal", "Obese", "…
## $ `Blood Pressure`          <chr> "126/83", "125/80", "125/80", "140/90", "140…
## $ `Heart Rate`              <dbl> 77, 75, 75, 85, 85, 85, 82, 70, 70, 70, 70, …
## $ `Daily Steps`             <dbl> 4200, 10000, 10000, 3000, 3000, 3000, 3500, …
## $ `Sleep Disorder`          <dbl> 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,…

skim(sueno)

Data summary

Name	sueno
Number of rows	374
Number of columns	13
_______________________
Column type frequency:
character	5
numeric	8
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
Gender	0	1	4	6	0	2	0
Occupation	0	1	5	20	0	11	0
Sleep Duration	0	1	1	3	0	27	0
BMI Category	0	1	5	13	0	4	0
Blood Pressure	0	1	6	6	0	25	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Person ID	0	1	187.50	108.11	1	94.25	187.5	280.75	374	▇▇▇▇▇
Age	0	1	42.18	8.67	27	35.25	43.0	50.00	59	▆▆▇▃▅
Quality of Sleep	0	1	7.31	1.20	4	6.00	7.0	8.00	9	▁▇▆▇▅
Physical Activity Level	0	1	59.17	20.83	30	45.00	60.0	75.00	90	▇▇▇▇▇
Stress Level	0	1	5.39	1.77	3	4.00	5.0	7.00	8	▇▃▂▃▃
Heart Rate	0	1	70.17	4.14	65	68.00	70.0	72.00	86	▇▇▂▁▁
Daily Steps	0	1	6816.84	1617.92	3000	5600.00	7000.0	8000.00	10000	▁▅▇▆▂
Sleep Disorder	0	1	0.41	0.49	0	0.00	0.0	1.00	1	▇▁▁▁▆

Interpretación:

• La base tiene 374 registros y 13 columnas, sin datos vacíos.
• Sleep Duration aparece como texto (chr) y Blood Pressure viene como “126/83”. Ambas deben convertirse a números.
• Sleep Disorder ya viene codificada como 0 = sin trastorno y 1 = con trastorno (insomnio o apnea).

# ¿Cuántos registros son copias exactas de otro (sin contar el ID)?
sum(duplicated(select(sueno, -`Person ID`)))

## [1] 248

Hallazgo de calidad de datos: 248 de los 374 registros son copias exactas de otro trabajador (mismos valores en todas las variables). Esto es típico de bases de práctica generadas artificialmente y se tendrá en cuenta al evaluar el modelo (sección 3.5.2).

3.2.2 Construcción, selección y limpieza de variables

base_modelo <- sueno %>%
  select(Gender, Age, `Sleep Duration`, `Quality of Sleep`, `Physical Activity Level`,
         `Stress Level`, `BMI Category`, `Blood Pressure`, `Heart Rate`, `Daily Steps`,
         `Sleep Disorder`) %>%
  mutate(
    sexo       = factor(Gender, levels = c("Female", "Male"), labels = c("Mujer", "Hombre")),
    edad       = as.numeric(Age),
    duracion   = as.numeric(`Sleep Duration`),            # viene como texto en el Excel
    calidad    = as.numeric(`Quality of Sleep`),
    actividad  = as.numeric(`Physical Activity Level`),    # minutos de actividad al día
    estres     = as.numeric(`Stress Level`),               # escala 1 a 10
    # "Normal" y "Normal Weight" son lo mismo; sobrepeso y obesidad se agrupan
    imc        = factor(ifelse(`BMI Category` %in% c("Normal", "Normal Weight"),
                               "Normal", "Sobrepeso/Obesidad")),
    sistolica  = as.numeric(sub("/.*", "", `Blood Pressure`)),   # número antes de "/"
    diastolica = as.numeric(sub(".*/", "", `Blood Pressure`)),   # número después de "/"
    fc         = as.numeric(`Heart Rate`),                 # frecuencia cardíaca
    pasos      = as.numeric(`Daily Steps`),
    trastorno  = factor(`Sleep Disorder`, levels = c(0, 1), labels = c("No", "Si"))
  ) %>%
  select(trastorno, calidad, sexo, edad, duracion, actividad, estres, imc,
         sistolica, diastolica, fc, pasos) %>%
  drop_na()

dim(base_modelo)

## [1] 374  12

summary(base_modelo)

##  trastorno    calidad          sexo          edad          duracion    
##  No:219    Min.   :4.000   Mujer :185   Min.   :27.00   Min.   :5.800  
##  Si:155    1st Qu.:6.000   Hombre:189   1st Qu.:35.25   1st Qu.:6.400  
##            Median :7.000                Median :43.00   Median :7.200  
##            Mean   :7.313                Mean   :42.18   Mean   :7.132  
##            3rd Qu.:8.000                3rd Qu.:50.00   3rd Qu.:7.800  
##            Max.   :9.000                Max.   :59.00   Max.   :8.500  
##    actividad         estres                      imc        sistolica    
##  Min.   :30.00   Min.   :3.000   Normal            :216   Min.   :115.0  
##  1st Qu.:45.00   1st Qu.:4.000   Sobrepeso/Obesidad:158   1st Qu.:125.0  
##  Median :60.00   Median :5.000                            Median :130.0  
##  Mean   :59.17   Mean   :5.385                            Mean   :128.6  
##  3rd Qu.:75.00   3rd Qu.:7.000                            3rd Qu.:135.0  
##  Max.   :90.00   Max.   :8.000                            Max.   :142.0  
##    diastolica          fc            pasos      
##  Min.   :75.00   Min.   :65.00   Min.   : 3000  
##  1st Qu.:80.00   1st Qu.:68.00   1st Qu.: 5600  
##  Median :85.00   Median :70.00   Median : 7000  
##  Mean   :84.65   Mean   :70.17   Mean   : 6817  
##  3rd Qu.:90.00   3rd Qu.:72.00   3rd Qu.: 8000  
##  Max.   :95.00   Max.   :86.00   Max.   :10000

Variables candidatas y justificación:

Variable	Tipo	Por qué se incluye
trastorno	Categórica Sí/No (respuesta del modelo principal)	Es lo que se quiere detectar.
calidad	Numérica 1–10 (respuesta del modelo complementario)	Indicador de bienestar que se quiere mejorar.
edad, sexo	Numérica / categórica	La apnea aumenta con la edad y es más frecuente en hombres.
duracion	Numérica (horas)	Dormir poco es síntoma típico del insomnio.
estres	Numérica 1–10	El estrés laboral se asocia con insomnio.
actividad, pasos	Numéricas	El sedentarismo empeora la calidad del sueño.
imc	Categórica	El sobrepeso y la obesidad son el principal factor de riesgo de apnea.
sistolica, fc	Numéricas	La hipertensión y la frecuencia cardíaca alta acompañan a la apnea.

Variables excluidas:

• Person ID: es solo un identificador.
• Occupation: tiene 11 categorías, varias con 1 a 4 personas, y el modelo debe servir para cualquier empresa.
• diastolica: va casi idéntica a la sistólica (correlación ≈ 0,97); usar ambas repetiría la misma información.

# Proporción de trabajadores con trastorno según IMC
base_modelo %>% count(imc, trastorno) %>% group_by(imc) %>%
  mutate(porcentaje = round(100 * n / sum(n), 1))

## # A tibble: 4 × 4
## # Groups:   imc [2]
##   imc                trastorno     n porcentaje
##   <fct>              <fct>     <int>      <dbl>
## 1 Normal             No          200       92.6
## 2 Normal             Si           16        7.4
## 3 Sobrepeso/Obesidad No           19       12  
## 4 Sobrepeso/Obesidad Si          139       88

# Promedio de las variables numéricas según presencia de trastorno
base_modelo %>% group_by(trastorno) %>%
  summarise(across(c(edad, duracion, estres, sistolica, fc, actividad), ~ round(mean(.x), 1)))

## # A tibble: 2 × 7
##   trastorno  edad duracion estres sistolica    fc actividad
##   <fct>     <dbl>    <dbl>  <dbl>     <dbl> <dbl>     <dbl>
## 1 No         39        7.4    5.1      124   69        57.9
## 2 Si         46.6      6.8    5.8      135.  71.8      60.9

Interpretación: el 88 % de los trabajadores con sobrepeso u obesidad tiene trastorno, contra el 7 % de los de peso normal. Quienes tienen trastorno también tienen mayor presión sistólica (≈ 135 vs. 124 mmHg), mayor edad y menos horas de sueño. Estas diferencias anticipan que el IMC y la presión arterial serán las variables con más peso.

3.3 Esquema conceptual del modelo

3.3.1 Diagrama

grViz("
digraph modelo {
  graph [rankdir = LR, fontsize = 11]
  node  [shape = box, style = 'rounded,filled', fillcolor = '#EAF2FB', fontname = Helvetica]

  subgraph cluster_x {
    label = 'Datos del examen ocupacional (X)'; style = dashed
    dem  [label = 'Edad y sexo']
    imc  [label = 'IMC']
    pa   [label = 'Presión sistólica\\ny frecuencia cardíaca']
    hab  [label = 'Horas de sueño, estrés,\\nactividad y pasos']
  }

  glm  [label = 'Regresión logística\\nglm(family = binomial)', shape = ellipse, fillcolor = '#D5F5E3']
  lm   [label = 'Regresión lineal\\nlm()', shape = ellipse, fillcolor = '#D5F5E3']
  y1   [label = 'Probabilidad de\\ntrastorno del sueño', fillcolor = '#FADBD8']
  y2   [label = 'Calidad del sueño\\n(1 a 10)', fillcolor = '#FADBD8']
  dec  [label = 'Priorizar remisión a\\npolisomnografía', shape = note, fillcolor = '#F9E79F']
  prog [label = 'Diseñar programa\\nde bienestar', shape = note, fillcolor = '#F9E79F']

  {dem imc pa hab} -> glm -> y1 -> dec
  {dem imc pa hab} -> lm  -> y2 -> prog
}
")

3.3.2 Pseudocódigo

1. Leer el Excel, convertir texto a números y separar la presión arterial.
2. Seleccionar variables y eliminar vacíos.
3. Dividir: 75 % entrenamiento / 25 % prueba (manteniendo la misma proporción de casos con trastorno).
4. Modelo principal: P(trastorno = Sí) = logística(b0 + b1*edad + b2*sexo + ... + b9*pasos)
5. Si P ≥ 0,5 → "remitir a estudio"; si no → "control rutinario".
6. Evaluar AUC, sensibilidad y especificidad en entrenamiento y prueba.
7. Modelo complementario: calidad = b0 + b1*duración + b2*estrés + ... + error
8. Diagnosticar residuos (error = real − predicho): ¿varianza constante? ¿forma normal?
9. Mejoras: evaluar sin registros duplicados y comparar con bosque aleatorio.

3.4 Implementación en RStudio

3.4.1 División en entrenamiento (75 %) y prueba (25 %)

division      <- initial_split(base_modelo, prop = 0.75, strata = trastorno)  # strata: misma proporción de "Sí" en ambos
entrenamiento <- training(division)
prueba        <- testing(division)

c(entrenamiento = nrow(entrenamiento), prueba = nrow(prueba))

## entrenamiento        prueba 
##           280            94

prop.table(table(entrenamiento$trastorno)); prop.table(table(prueba$trastorno))

## 
##        No        Si 
## 0.5857143 0.4142857

## 
##        No        Si 
## 0.5851064 0.4148936

Interpretación: el modelo aprende con 280 trabajadores y se califica con 94 que no vio. En ambos grupos, cerca del 41 % tiene trastorno, por lo que la comparación es justa.

3.4.2 Modelo principal: regresión logística

formula_clas <- trastorno ~ edad + sexo + duracion + estres + actividad + imc + sistolica + fc + pasos

modelo_log <- glm(formula_clas, data = entrenamiento, family = binomial)
# Equivalente con caret: caret::train(formula_clas, data = entrenamiento, method = "glm", family = "binomial")

tidy(modelo_log) %>%
  mutate(odds_ratio = exp(estimate)) %>%            # cuántas veces cambia el riesgo
  mutate(across(where(is.numeric), ~ round(.x, 3)))

## # A tibble: 10 × 6
##    term                  estimate std.error statistic p.value odds_ratio
##    <chr>                    <dbl>     <dbl>     <dbl>   <dbl>      <dbl>
##  1 (Intercept)            -68.9      18.1      -3.82    0          0    
##  2 edad                    -0.054     0.08     -0.673   0.501      0.948
##  3 sexoHombre              -0.398     0.696    -0.572   0.567      0.671
##  4 duracion                 1.95      1.24      1.57    0.117      7.04 
##  5 estres                   1.23      0.506     2.43    0.015      3.42 
##  6 actividad               -0.041     0.032    -1.27    0.203      0.96 
##  7 imcSobrepeso/Obesidad    2.66      0.907     2.93    0.003     14.3  
##  8 sistolica                0.43      0.106     4.06    0          1.54 
##  9 fc                      -0.062     0.125    -0.494   0.621      0.94 
## 10 pasos                    0         0         0.019   0.984      1

Interpretación (odds_ratio): razón de momios, es decir, cuántas veces se multiplica el riesgo de tener trastorno cuando la variable aumenta en 1 unidad, manteniendo lo demás igual. Solo se interpretan las variables con valor p < 0,05.

• IMC sobrepeso/obesidad (OR ≈ 14): el riesgo es unas 14 veces mayor que con peso normal. Es el factor más fuerte.
• Presión sistólica (OR ≈ 1,5): cada mmHg adicional aumenta el riesgo cerca de un 50 %.
• Estrés (OR ≈ 3,4): cada punto adicional en la escala de estrés multiplica el riesgo por más de 3.
• Edad, sexo, horas de sueño, actividad, frecuencia cardíaca y pasos no son significativos cuando ya se conocen IMC, presión y estrés.

3.4.3 Desempeño del modelo principal: entrenamiento vs. prueba

umbral <- 0.5
prob_ent <- predict(modelo_log, newdata = entrenamiento, type = "response")  # probabilidad de "Si"
prob_pru <- predict(modelo_log, newdata = prueba,        type = "response")

metricas_clas <- function(real, prob) {
  pred <- factor(ifelse(prob >= umbral, "Si", "No"), levels = c("No", "Si"))
  c(AUC           = roc_auc_vec(real, prob, event_level = "second"),
    Sensibilidad  = sens_vec(real, pred, event_level = "second"),
    Especificidad = spec_vec(real, pred, event_level = "second"),
    Exactitud     = accuracy_vec(real, pred))
}

tabla_log <- rbind(Entrenamiento = metricas_clas(entrenamiento$trastorno, prob_ent),
                   Prueba        = metricas_clas(prueba$trastorno, prob_pru))
round(tabla_log, 3)

##                 AUC Sensibilidad Especificidad Exactitud
## Entrenamiento 0.949        0.922         0.963     0.946
## Prueba        0.921        0.923         0.945     0.936

# Matriz de confusión en prueba (filas = real, columnas = predicho)
table(Real = prueba$trastorno,
      Predicho = factor(ifelse(prob_pru >= umbral, "Si", "No"), levels = c("No", "Si")))

##     Predicho
## Real No Si
##   No 52  3
##   Si  3 36

Interpretación:

• AUC en prueba ≈ 0,92, por encima de la meta de 0,85. La brecha con entrenamiento (≈ 0,95) es pequeña.
• Sensibilidad ≈ 0,92: de cada 10 trabajadores con trastorno, el modelo detecta 9.
• Especificidad ≈ 0,95: casi no envía a estudio a personas sanas.
• En la matriz de confusión, los errores son pocos: unos 3 casos no detectados y 3 remisiones innecesarias entre 94 trabajadores.

# Curva ROC: muestra sensibilidad y falsos positivos para cada umbral posible
umbrales <- seq(0, 1, by = 0.01)
curva_roc <- data.frame(
  umbral        = umbrales,
  sensibilidad  = sapply(umbrales, function(u) mean(prob_pru[prueba$trastorno == "Si"] >= u)),
  falsos_pos    = sapply(umbrales, function(u) mean(prob_pru[prueba$trastorno == "No"] >= u))
)

g_roc <- ggplot(curva_roc, aes(x = falsos_pos, y = sensibilidad, text = paste("Umbral:", umbral))) +
  geom_path(aes(group = 1), color = "steelblue", linewidth = 1) +
  geom_abline(linetype = "dashed", color = "grey50") +
  labs(title = "Curva ROC del modelo logístico (conjunto de prueba)",
       x = "1 − Especificidad (proporción de falsos positivos)", y = "Sensibilidad") +
  theme_minimal()
ggplotly(g_roc, tooltip = c("text", "x", "y"))

Comentario: la curva sube pegada al eje izquierdo y se aleja de la diagonal (que representa adivinar al azar). Al pasar el mouse se ve el umbral de cada punto. Si la IPS quiere detectar aún más casos, puede bajar el umbral (por ejemplo, a 0,3), aceptando más remisiones innecesarias.

3.4.4 Visualización 3D interactiva

plot_ly(base_modelo,
        x = ~jitter(sistolica, 1), y = ~jitter(estres, 0.6), z = ~jitter(edad, 1),   # jitter: separa puntos repetidos
        color = ~trastorno, colors = c("#2E86C1", "#E74C3C"),
        symbol = ~imc, symbols = c("circle", "diamond"),
        type = "scatter3d", mode = "markers", marker = list(size = 4, opacity = 0.7)) %>%
  layout(title = "Trabajadores según presión sistólica, estrés y edad",
         scene = list(xaxis = list(title = "Presión sistólica (mmHg)"),
                      yaxis = list(title = "Estrés (1-10)"),
                      zaxis = list(title = "Edad")))

Comentario: los trabajadores con trastorno (rojo) se concentran en la zona de presión sistólica alta, y la mayoría tiene forma de rombo (sobrepeso u obesidad). Los de azul se agrupan con presión por debajo de 130 mmHg.

# Superficie de probabilidad predicha según presión sistólica y estrés,
# para un trabajador con peso normal y otro con sobrepeso/obesidad (demás variables en la mediana)
rejilla <- expand.grid(sistolica = seq(115, 142, length.out = 30),
                       estres    = seq(3, 8, length.out = 30))
valores_tipicos <- base_modelo %>%
  summarise(across(c(edad, duracion, actividad, fc, pasos), median))
rejilla <- cbind(rejilla, valores_tipicos, sexo = factor("Mujer", levels = levels(base_modelo$sexo)))

prob_normal    <- predict(modelo_log, type = "response",
                          newdata = mutate(rejilla, imc = factor("Normal", levels = levels(base_modelo$imc))))
prob_sobrepeso <- predict(modelo_log, type = "response",
                          newdata = mutate(rejilla, imc = factor("Sobrepeso/Obesidad", levels = levels(base_modelo$imc))))

x_sis <- unique(rejilla$sistolica); y_est <- unique(rejilla$estres)
plot_ly() %>%
  add_surface(x = x_sis, y = y_est, z = t(matrix(prob_normal, nrow = 30)),
              colorscale = list(c(0, 1), c("#D6EAF8", "#1B4F72")), showscale = FALSE,
              name = "IMC normal", opacity = 0.9) %>%
  add_surface(x = x_sis, y = y_est, z = t(matrix(prob_sobrepeso, nrow = 30)),
              colorscale = list(c(0, 1), c("#FADBD8", "#922B21")), showscale = FALSE,
              name = "Sobrepeso/obesidad", opacity = 0.9) %>%
  layout(title = "Probabilidad de trastorno (azul: IMC normal; rojo: sobrepeso)",
         scene = list(xaxis = list(title = "Presión sistólica"),
                      yaxis = list(title = "Estrés"),
                      zaxis = list(title = "Probabilidad", range = c(0, 1)),
                      camera = list(eye = list(x = 1.6, y = -1.6, z = 0.9))))

Comentario: la superficie roja (sobrepeso/obesidad) sube mucho antes que la azul. Con sobrepeso, una presión de 130 mmHg ya implica alta probabilidad de trastorno; con peso normal, la probabilidad solo aumenta cuando se combinan presión alta y estrés alto. Esta gráfica se puede usar como regla práctica de remisión en el examen ocupacional.

3.4.5 Modelo complementario: regresión lineal de la calidad del sueño

formula_reg <- calidad ~ duracion + estres + edad + fc + actividad + pasos + sexo + imc

modelo_lin <- lm(formula_reg, data = entrenamiento)
tidy(modelo_lin) %>% mutate(across(where(is.numeric), ~ round(.x, 4)))

## # A tibble: 9 × 5
##   term                  estimate std.error statistic p.value
##   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)             7.46      0.639      11.7   0     
## 2 duracion                0.141     0.0519      2.72  0.007 
## 3 estres                 -0.417     0.0255    -16.3   0     
## 4 edad                    0.0534    0.0037     14.6   0     
## 5 fc                     -0.0215    0.0078     -2.76  0.0061
## 6 actividad               0.0044    0.0019      2.39  0.0177
## 7 pasos                   0         0           1.62  0.107 
## 8 sexoHombre              0.244     0.0492      4.96  0     
## 9 imcSobrepeso/Obesidad  -0.750     0.0613    -12.2   0

glance(modelo_lin) %>% select(r.squared, adj.r.squared, sigma) %>% round(3)

## # A tibble: 1 × 3
##   r.squared adj.r.squared sigma
##       <dbl>         <dbl> <dbl>
## 1     0.949         0.947 0.279

pred_ent_lin <- predict(modelo_lin, newdata = entrenamiento)
pred_pru_lin <- predict(modelo_lin, newdata = prueba)

tabla_lin <- rbind(
  Entrenamiento = c(RMSE = rmse_vec(entrenamiento$calidad, pred_ent_lin),
                    R2   = rsq_vec(entrenamiento$calidad, pred_ent_lin)),
  Prueba        = c(RMSE = rmse_vec(prueba$calidad, pred_pru_lin),
                    R2   = rsq_vec(prueba$calidad, pred_pru_lin)))
round(tabla_lin, 3)

##                RMSE    R2
## Entrenamiento 0.275 0.949
## Prueba        0.344 0.912

Interpretación:

• Estrés: cada punto adicional reduce la calidad del sueño en ≈ 0,42 puntos. Es la palanca modificable más importante.
• Sobrepeso/obesidad: reduce la calidad en ≈ 0,75 puntos.
• Actividad física: tiene un efecto positivo pero pequeño.
• Ajuste: el modelo explica cerca del 95 % de la variación (R²). El RMSE en prueba ≈ 0,34 puntos cumple la meta (≤ 0,5).

3.4.6 Diagnóstico preliminar de supuestos (modelo lineal)

La regresión lineal supone que los residuos (error = calidad real − calidad predicha):

1. tienen varianza constante (homocedasticidad: el error es del mismo tamaño para calidades bajas y altas), y
2. siguen una forma aproximadamente normal (campana centrada en 0).

broom::augment() agrega a cada fila su valor predicho (.fitted) y su residuo (.resid). ggplot2 no tiene una función geom_residuals(), por lo que los residuos se grafican con geom_point(), geom_histogram() y stat_qq().

diagnostico <- augment(modelo_lin)

g_res <- ggplot(diagnostico, aes(x = .fitted, y = .resid)) +
  geom_jitter(alpha = 0.5, color = "steelblue", width = 0.03, height = 0.03) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  geom_smooth(se = FALSE, color = "darkorange", method = "loess") +
  labs(title = "Residuos vs. valores predichos (calidad del sueño)",
       x = "Calidad predicha", y = "Residuo (puntos)") +
  theme_minimal()
ggplotly(g_res)

g_hist <- ggplot(diagnostico, aes(x = .resid)) +
  geom_histogram(bins = 25, fill = "steelblue", color = "white") +
  labs(title = "Distribución de los residuos", x = "Residuo (puntos)", y = "Frecuencia") +
  theme_minimal()
ggplotly(g_hist)

ggplot(diagnostico, aes(sample = .resid)) +
  stat_qq(alpha = 0.5) + stat_qq_line(color = "red") +
  labs(title = "Gráfico Q-Q de los residuos",
       x = "Cuantiles teóricos (normal)", y = "Cuantiles de los residuos") +
  theme_minimal()

shapiro.test(diagnostico$.resid)   # prueba de normalidad: p < 0,05 indica que no es normal

## 
##  Shapiro-Wilk normality test
## 
## data:  diagnostico$.resid
## W = 0.96354, p-value = 1.635e-06

Interpretación:

• Homocedasticidad: los residuos están centrados en 0, pero forman franjas, porque la calidad del sueño solo toma valores enteros (4, 5, 6…). La dispersión es algo mayor en las calidades bajas. La varianza constante se cumple solo de forma aproximada.
• Normalidad: el histograma se parece a una campana, pero en el gráfico Q-Q los extremos se separan de la línea y la prueba de Shapiro-Wilk da p < 0,05. Los residuos no son exactamente normales.
• Conclusión del diagnóstico: el modelo sirve para describir qué factores influyen. Como la respuesta es una escala de enteros, un modelo para datos ordinales (categorías ordenadas) sería técnicamente más apropiado; ese es tema de un curso posterior.

3.5 Interpretación de resultados

3.5.1 Resumen de métricas frente a las metas

data.frame(
  Metrica = c("AUC (prueba)", "Sensibilidad (prueba)", "Especificidad (prueba)",
              "Brecha AUC entrenamiento - prueba", "RMSE calidad (prueba)"),
  Resultado = round(c(tabla_log["Prueba", "AUC"], tabla_log["Prueba", "Sensibilidad"],
                      tabla_log["Prueba", "Especificidad"],
                      tabla_log["Entrenamiento", "AUC"] - tabla_log["Prueba", "AUC"],
                      tabla_lin["Prueba", "RMSE"]), 3),
  Meta = c(">= 0,85", ">= 0,80", ">= 0,75", "<= 0,05", "<= 0,50")
)

##                             Metrica Resultado    Meta
## 1                      AUC (prueba)     0.921 >= 0,85
## 2             Sensibilidad (prueba)     0.923 >= 0,80
## 3            Especificidad (prueba)     0.945 >= 0,75
## 4 Brecha AUC entrenamiento - prueba     0.028 <= 0,05
## 5             RMSE calidad (prueba)     0.344 <= 0,50

Interpretación: con la base completa, los dos modelos cumplen todas las metas. Antes de recomendar su uso, hay que revisar si el resultado es demasiado optimista por los registros duplicados.

3.5.2 Acción de mejora 1 (implementada): evaluación sin registros duplicados

# ¿Qué porcentaje de los trabajadores de prueba tiene una "copia exacta" en entrenamiento?
clave <- function(df) do.call(paste, df)
round(100 * mean(clave(prueba) %in% clave(entrenamiento)), 1)

## [1] 75.5

Interpretación: cerca del 75 % de los trabajadores de prueba (3 de cada 4) ya aparecen idénticos en entrenamiento. El modelo está siendo calificado en gran parte con “preguntas que ya vio”, lo que se conoce como fuga de información. Para medir el desempeño real, se repite todo con perfiles únicos.

base_unica <- distinct(base_modelo)     # deja un solo registro por perfil
nrow(base_unica)

## [1] 122

division_u  <- initial_split(base_unica, prop = 0.75, strata = trastorno)
ent_u <- training(division_u); pru_u <- testing(division_u)

modelo_log_u <- glm(formula_clas, data = ent_u, family = binomial)

tabla_unica <- rbind(
  Entrenamiento = metricas_clas(ent_u$trastorno, predict(modelo_log_u, ent_u, type = "response")),
  Prueba        = metricas_clas(pru_u$trastorno, predict(modelo_log_u, pru_u, type = "response")))
round(tabla_unica, 3)

##                 AUC Sensibilidad Especificidad Exactitud
## Entrenamiento 0.894        0.816         0.849     0.835
## Prueba        0.825        0.846         0.778     0.806

Interpretación:

• Solo quedan 122 perfiles únicos, así que el conjunto de prueba queda de 31 trabajadores y las métricas varían más según cuáles caigan en él.
• Sin duplicados, el AUC en prueba baja de 0,92 a ≈ 0,83, la sensibilidad a ≈ 0,85 y la especificidad a ≈ 0,78. Esta es la estimación más realista del desempeño en trabajadores nuevos.
• Con esta evaluación, el modelo queda levemente por debajo de la meta de AUC (0,85), aunque cumple sensibilidad y especificidad. Sirve como herramienta de priorización, pero no reemplaza la valoración clínica y necesita más datos reales (sección 3.5.4).

3.5.3 Acción de mejora 2 (implementada): probar otro algoritmo

Un bosque aleatorio combina cientos de árboles de decisión y puede capturar relaciones no lineales. Se entrena con las mismas variables y se califica con el mismo conjunto de prueba de la sección 3.4.

modelo_rf <- randomForest(formula_clas, data = entrenamiento, ntree = 500, importance = TRUE)
# Equivalente con caret: caret::train(formula_clas, data = entrenamiento, method = "rf")

prob_ent_rf <- predict(modelo_rf, newdata = entrenamiento, type = "prob")[, "Si"]
prob_pru_rf <- predict(modelo_rf, newdata = prueba,        type = "prob")[, "Si"]

comparacion <- data.frame(
  Modelo   = rep(c("Regresión logística", "Bosque aleatorio"), each = 2),
  Conjunto = rep(c("Entrenamiento", "Prueba"), 2),
  rbind(metricas_clas(entrenamiento$trastorno, prob_ent),
        metricas_clas(prueba$trastorno, prob_pru),
        metricas_clas(entrenamiento$trastorno, prob_ent_rf),
        metricas_clas(prueba$trastorno, prob_pru_rf)))
comparacion %>% mutate(across(where(is.numeric), ~ round(.x, 3)))

##                Modelo      Conjunto   AUC Sensibilidad Especificidad Exactitud
## 1 Regresión logística Entrenamiento 0.949        0.922         0.963     0.946
## 2 Regresión logística        Prueba 0.921        0.923         0.945     0.936
## 3    Bosque aleatorio Entrenamiento 0.980        0.922         0.963     0.946
## 4    Bosque aleatorio        Prueba 0.916        0.923         0.945     0.936

importancia <- data.frame(variable = rownames(importance(modelo_rf)),
                          importancia = importance(modelo_rf)[, "MeanDecreaseGini"]) %>%
  arrange(importancia)

g_imp <- ggplot(importancia, aes(x = reorder(variable, importancia), y = importancia)) +
  geom_col(fill = "darkgreen") + coord_flip() +
  labs(title = "Importancia de las variables (bosque aleatorio)",
       x = "", y = "Disminución media de impureza (Gini)") +
  theme_minimal()
ggplotly(g_imp)

Interpretación:

• El bosque aleatorio tiene un AUC muy alto en entrenamiento, pero en prueba no supera a la regresión logística. La brecha más grande indica que tiende a “memorizar” los datos.
• Con pocos datos, el modelo más simple (logístico) es preferible, porque además es más fácil de explicar al equipo médico y a la gerencia.
• El gráfico de importancia confirma que IMC y presión sistólica son las variables que más aportan, seguidas de la edad y la duración del sueño.

3.5.4 Otras acciones de mejora propuestas

1. Recolección de datos adicionales:
   • Aumentar el número de trabajadores reales y no repetidos.
   • Volver a separar el tipo de trastorno (insomnio vs. apnea), porque tienen causas y tratamientos distintos.
   • Agregar variables de tamizaje validadas para apnea, como ronquido, perímetro del cuello y el cuestionario STOP-Bang.
   • Registrar si el trabajador hace turnos nocturnos.
2. Ingeniería de variables:
   • Clasificar la presión arterial según las guías clínicas (normal, elevada, hipertensión grado 1 y 2).
   • Agrupar las ocupaciones por tipo de jornada (asistencial con turnos, administrativa, comercial) para incluirlas en el modelo sin categorías tan pequeñas.
3. Ajuste del umbral de decisión: en lugar de 0,5, elegir el umbral con la curva ROC según los cupos de polisomnografía disponibles. Un umbral más bajo detecta más casos, pero genera más remisiones.

3.5.5 Conclusiones

1. Con datos de un examen ocupacional básico, la regresión logística identifica a los trabajadores con trastorno del sueño con un AUC ≈ 0,92 en la base completa y ≈ 0,83 con perfiles únicos (evaluación más realista).
2. El sobrepeso/obesidad, la presión sistólica y el estrés son los factores que más aumentan el riesgo. Los programas de prevención de la IPS deben priorizar el control de peso, de presión arterial y del estrés laboral.
3. El estrés también es el factor que más reduce la calidad del sueño, lo que refuerza la importancia de intervenir el riesgo psicosocial.
4. Los registros duplicados inflan las métricas. La evaluación con perfiles únicos es la referencia más honesta, y el modelo debe validarse con datos reales de la población trabajadora antes de usarse para decidir remisiones.