R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

# ============================================================
# I. EXPLORACIÓN DE LA BASE DE DATOS 
# ============================================================

datos <- read_excel("salud_cardiovascular.xlsx")
datos <- as.data.frame(datos)
dim(datos) 

[1] 500  10

names(datos)

 [1] "id"                    "edad"                  "sexo"                 
 [4] "imc"                   "presion_sistolica"     "glucosa_ayuno"        
 [7] "actividad_fisica"      "fumador"               "riesgo_cardiovascular"
[10] "riesgo_alto"          

head(datos)

# Revisamos la estructura de la base como qué variables son numéricas, cuáles son texto, factores, etc. 
str(datos)

'data.frame':   500 obs. of  10 variables:
 $ id                   : num  1 2 3 4 5 6 7 8 9 10 ...
 $ edad                 : num  19 24 44 74 31 31 47 70 46 69 ...
 $ sexo                 : chr  "F" "M" "F" "F" ...
 $ imc                  : num  23.3 17.8 25 33.1 28.9 33.1 25.2 31.2 26.4 25.4 ...
 $ presion_sistolica    : num  102 138 116 149 80 118 132 123 117 118 ...
 $ glucosa_ayuno        : num  73 98 84 91 84 97 104 115 92 76 ...
 $ actividad_fisica     : chr  "activo" "moderado" "activo" "moderado" ...
 $ fumador              : chr  "si" "no" "no" "no" ...
 $ riesgo_cardiovascular: num  18.4 20.6 17.1 73.8 28.6 27.1 42 63.7 25.7 73.4 ...
 $ riesgo_alto          : chr  "bajo" "bajo" "bajo" "alto" ...

# Resumen general de todas las variables.
summary(datos)

       id             edad          sexo                imc       
 Min.   :  1.0   Min.   :18.0   Length:500         Min.   :16.00  
 1st Qu.:125.8   1st Qu.:33.0   Class :character   1st Qu.:23.90  
 Median :250.5   Median :48.0   Mode  :character   Median :27.15  
 Mean   :250.5   Mean   :48.4                      Mean   :27.26  
 3rd Qu.:375.2   3rd Qu.:64.0                      3rd Qu.:30.70  
 Max.   :500.0   Max.   :80.0                      Max.   :41.40  
 presion_sistolica glucosa_ayuno   actividad_fisica     fumador         
 Min.   : 80.0     Min.   : 60.0   Length:500         Length:500        
 1st Qu.:109.0     1st Qu.: 81.0   Class :character   Class :character  
 Median :120.0     Median : 95.0   Mode  :character   Mode  :character  
 Mean   :120.0     Mean   : 95.2                                        
 3rd Qu.:130.2     3rd Qu.:110.0                                        
 Max.   :159.0     Max.   :157.0                                        
 riesgo_cardiovascular riesgo_alto       
 Min.   :  0.00        Length:500        
 1st Qu.: 26.00        Class :character  
 Median : 37.75        Mode  :character  
 Mean   : 40.12                          
 3rd Qu.: 53.23                          
 Max.   :100.00                          

# Para ver cuántos datos faltantes hay por variable.
colSums(is.na(datos))

                   id                  edad                  sexo 
                    0                     0                     0 
                  imc     presion_sistolica         glucosa_ayuno 
                    0                     0                     0 
     actividad_fisica               fumador riesgo_cardiovascular 
                    0                     0                     0 
          riesgo_alto 
                    0 

# ============================================================
# II. DEFINIR VARIABLE DEPENDIENTE E INDEPENDIENTES
# ============================================================

# Para hacerla factor para que R la trate como categoría.
datos$riesgo_alto <- factor(datos$riesgo_alto)

# También las hago a factor las variables categóricas independientes
datos$sexo <- factor(datos$sexo)
datos$actividad_fisica <- factor(datos$actividad_fisica)
datos$fumador <- factor(datos$fumador)

# cuántas personas hay en cada categoría de riesgo.
table(datos$riesgo_alto)

alto bajo 
 149  351 

# el porcentaje de personas en cada categoría.
prop.table(table(datos$riesgo_alto))

 alto  bajo 
0.298 0.702 

# ============================================================
# III. VARIABLES QUE SÍ SE USARÁN COMO INDEPENDIENTES
# ============================================================

# Variables numéricas de salud.
variables_numericas <- c(
  "edad",
  "imc",
  "presion_sistolica",
  "glucosa_ayuno"
)

# Variables categóricas.
variables_categoricas <- c(
  "sexo",
  "actividad_fisica",
  "fumador"
)

# Solo para ver que quedaron bien guardadas. 
variables_numericas

[1] "edad"              "imc"               "presion_sistolica"
[4] "glucosa_ayuno"    

variables_categoricas

[1] "sexo"             "actividad_fisica" "fumador"         

# ============================================================
# IV. ANÁLISIS VISUAL 1: DISTRIBUCIÓN DE LA VARIABLE DEPENDIENTE
# ============================================================
barplot(
  table(datos$riesgo_alto),
  main = "Distribución del riesgo cardiovascular",
  xlab = "Nivel de riesgo",
  ylab = "Número de personas",
  col = "purple"
)

# ============================================================
# V. ANÁLISIS VISUAL 2: VARIABLES NUMÉRICAS VS RIESGO ALTO
# ============================================================

# Edad según nivel de riesgo.
boxplot(
  edad ~ riesgo_alto,
  data = datos,
  main = "Edad según riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Edad",
  col = "cyan"
)

# IMC según nivel de riesgo.
boxplot(
  imc ~ riesgo_alto,
  data = datos,
  main = "IMC según riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "IMC",
  col = "magenta"
)

# Presión sistólica según nivel de riesgo.
boxplot(
  presion_sistolica ~ riesgo_alto,
  data = datos,
  main = "Presión sistólica según riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Presión sistólica",
  col = "yellow"
)

# Glucosa en ayuno según nivel de riesgo.
boxplot(
  glucosa_ayuno ~ riesgo_alto,
  data = datos,
  main = "Glucosa en ayuno según riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Glucosa en ayuno",
  col = "orange"
)

# ============================================================
# VI. ANÁLISIS VISUAL 3: VARIABLES CATEGÓRICAS VS RIESGO ALTO
# ============================================================

# *****************************
# Sexo y riesgo cardiovascular
# *****************************

# Tabla de conteos.
table(datos$sexo, datos$riesgo_alto)

   
    alto bajo
  F   70  192
  M   79  159

# Tabla de porcentajes por sexo.
prop.table(table(datos$sexo, datos$riesgo_alto), 1)

   
         alto      bajo
  F 0.2671756 0.7328244
  M 0.3319328 0.6680672

# Gráfica de barras.
barplot(
  table(datos$sexo, datos$riesgo_alto),
  beside = TRUE,
  legend.text = TRUE,
  main = "Sexo y riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Número de personas",
  col = c("lightgreen", "darkgreen")
)

# ****************************************
# Actividad física y riesgo cardiovascular
# ****************************************

# Tabla de conteos.
table(datos$actividad_fisica, datos$riesgo_alto)

            
             alto bajo
  activo       23   73
  moderado     57  128
  sedentario   69  150

# Tabla de porcentajes por actividad física.
prop.table(table(datos$actividad_fisica, datos$riesgo_alto), 1)

            
                  alto      bajo
  activo     0.2395833 0.7604167
  moderado   0.3081081 0.6918919
  sedentario 0.3150685 0.6849315

# Gráfica de barras.
barplot(
  table(datos$actividad_fisica, datos$riesgo_alto),
  beside = TRUE,
  legend.text = TRUE,
  main = "Actividad física y riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Número de personas",
  col = c("skyblue", "steelblue", "darkblue")
)

# ********************************
# Fumador y riesgo cardiovascular
# ********************************

# Tabla de conteos.
table(datos$fumador, datos$riesgo_alto)

    
     alto bajo
  no  105  262
  si   44   89

# Tabla de porcentajes por fumador.
prop.table(table(datos$fumador, datos$riesgo_alto), 1)

    
          alto      bajo
  no 0.2861035 0.7138965
  si 0.3308271 0.6691729

# Gráfica de barras.
barplot(
  table(datos$fumador, datos$riesgo_alto),
  beside = TRUE, #para comparar mejor las categorías y que salgan pegaditas 
  legend.text = TRUE,
  main = "Fumador y riesgo cardiovascular",
  xlab = "Riesgo cardiovascular",
  ylab = "Número de personas",
  col = c("firebrick", "darkred")
)

# ================================================================================
# VII. EXPLORAR LA RELACIÓN ENTRE CADA VARIABLE INDEPENDIENTE Y LA VARIABLE DEPENDIENTE
# ================================================================================

# Creamos una versión numérica de la variable dependiente.
# riesgo_num vale 1 si la persona tiene riesgo alto.
# riesgo_num vale 0 si la persona tiene riesgo bajo.
datos$riesgo_num <- 0
datos$riesgo_num[datos$riesgo_alto == "alto"] <- 1

# Solo para revisar que sí se hizo la conversión 
table(datos$riesgo_alto, datos$riesgo_num) 

      
         0   1
  alto   0 149
  bajo 351   0

# **********************************************
# a. MODELO 1: RELACIÓN ENTRE EDAD Y RIESGO ALTO
# **********************************************

# Queremos ver cómo cambia el riesgo alto conforme cambia la edad.
modelo_loess_edad <- loess(riesgo_num ~ edad, data = datos, span = 0.75)

# Creamos una secuencia de edades para poder dibujar la curva suavizada.
x_edad <- seq(
  min(datos$edad),
  max(datos$edad),
  length.out = 100
) 

# Calculamos las predicciones del modelo LOESS para esas edades.
pred_edad <- predict(
  modelo_loess_edad,
  newdata = data.frame(edad = x_edad)
)

# Graficamos los puntos reales.
# Como riesgo_num solo vale 0 o 1, los puntos se verán abajo o arriba.
plot(
  datos$edad,
  datos$riesgo_num,
  main = "LOESS: Edad y riesgo cardiovascular alto",
  xlab = "Edad",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

# Agregamos la curva LOESS.
lines(
  x_edad,
  pred_edad,
  col = "cyan",
  lwd = 3
)

# ***************************************
# b. MODELO 2: IMC Y RIESGO ALTO CON LOESS
# ***************************************

# Ajustamos el modelo LOESS para IMC.
modelo_loess_imc <- loess(riesgo_num ~ imc, data = datos, span = 0.75)

# Creamos una secuencia de valores de IMC.
x_imc <- seq(
  min(datos$imc),
  max(datos$imc),
  length.out = 100
)

# Calculamos predicciones.
pred_imc <- predict(
  modelo_loess_imc,
  newdata = data.frame(imc = x_imc)
)

# Graficamos puntos y curva.
plot(
  datos$imc,
  datos$riesgo_num,
  main = "LOESS: IMC y riesgo cardiovascular alto",
  xlab = "IMC",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(
  x_imc,
  pred_imc,
  col = "magenta",
  lwd = 3
)

# *****************************************************
# c. MODELO 3: PRESIÓN SISTÓLICA Y RIESGO ALTO CON LOESS
# *****************************************************

# Ajustamos el modelo LOESS para presión sistólica.
modelo_loess_presion <- loess(riesgo_num ~ presion_sistolica, data = datos, span = 0.75)

# Creamos una secuencia de valores de presión.
x_presion <- seq(
  min(datos$presion_sistolica),
  max(datos$presion_sistolica),
  length.out = 100
)

# Calculamos predicciones.
pred_presion <- predict(
  modelo_loess_presion,
  newdata = data.frame(presion_sistolica = x_presion)
)

# Graficamos puntos y curva.
plot(
  datos$presion_sistolica,
  datos$riesgo_num,
  main = "LOESS: Presión sistólica y riesgo cardiovascular alto",
  xlab = "Presión sistólica",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(
  x_presion,
  pred_presion,
  col = "yellow",
  lwd = 3
)

# *****************************************************
# d. MODELO 4: GLUCOSA EN AYUNO Y RIESGO ALTO CON LOESS
# *****************************************************

# Ajustamos el modelo LOESS para glucosa en ayuno.
modelo_loess_glucosa <- loess(riesgo_num ~ glucosa_ayuno, data = datos, span = 0.75)

# Creamos una secuencia de valores de glucosa.
x_glucosa <- seq(
  min(datos$glucosa_ayuno),
  max(datos$glucosa_ayuno),
  length.out = 100
)

# Calculamos predicciones.
pred_glucosa <- predict(
  modelo_loess_glucosa,
  newdata = data.frame(glucosa_ayuno = x_glucosa)
)

# Graficamos puntos y curva.
plot(
  datos$glucosa_ayuno,
  datos$riesgo_num,
  main = "LOESS: Glucosa en ayuno y riesgo cardiovascular alto",
  xlab = "Glucosa en ayuno",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(
  x_glucosa,
  pred_glucosa,
  col = "orange",
  lwd = 3
)

# ============================================================
# VIII. REVISIÓN DE LA SENSIBILIDAD DEL SPAN POR VARIABLE
# ============================================================

# **********************************************
# a. MODELO 1: RELACIÓN ENTRE EDAD Y RIESGO ALTO
# **********************************************

# 0.4 = curva más flexible.
# 0.75 = curva intermedia, la que ya usamos.
# 0.9 = curva más suave.
loess_edad_04 <- loess(riesgo_num ~ edad, data = datos, span = 0.4)
loess_edad_075 <- loess(riesgo_num ~ edad, data = datos, span = 0.75)
loess_edad_09 <- loess(riesgo_num ~ edad, data = datos, span = 0.9)

# Secuencia de edades para graficar las curvas.
x_edad <- seq(min(datos$edad), max(datos$edad), length.out = 100)

# Predicciones para cada modelo.
pred_edad_04 <- predict(loess_edad_04, newdata = data.frame(edad = x_edad))
pred_edad_075 <- predict(loess_edad_075, newdata = data.frame(edad = x_edad))
pred_edad_09 <- predict(loess_edad_09, newdata = data.frame(edad = x_edad))

# Graficamos los puntos reales.
plot(
  datos$edad,
  datos$riesgo_num,
  main = "Comparación de span en LOESS: Edad",
  xlab = "Edad",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

# Agregamos las tres curvas.
lines(x_edad, pred_edad_04, col = "navy", lwd = 2)

lines(x_edad, pred_edad_075, col = "cyan", lwd = 3)
lines(x_edad, pred_edad_09, col = "brown", lwd = 2)

# Agregamos leyenda.
legend(
  "topleft",
  legend = c("span = 0.4", "span = 0.75", "span = 0.9"),
  col = c("navy", "cyan", "brown"),
  lwd = c(2, 3, 2)
)

# ***************************************
# b. MODELO 2: IMC Y RIESGO ALTO CON LOESS
# ***************************************

loess_imc_04 <- loess(riesgo_num ~ imc, data = datos, span = 0.4)
loess_imc_075 <- loess(riesgo_num ~ imc, data = datos, span = 0.75)
loess_imc_09 <- loess(riesgo_num ~ imc, data = datos, span = 0.9)

x_imc <- seq(min(datos$imc), max(datos$imc), length.out = 100)

pred_imc_04 <- predict(loess_imc_04, newdata = data.frame(imc = x_imc))
pred_imc_075 <- predict(loess_imc_075, newdata = data.frame(imc = x_imc))
pred_imc_09 <- predict(loess_imc_09, newdata = data.frame(imc = x_imc))

plot(
  datos$imc,
  datos$riesgo_num,
  main = "Comparación de span en LOESS: IMC",
  xlab = "IMC",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(x_imc, pred_imc_04, col = "navy", lwd = 2)

lines(x_imc, pred_imc_075, col = "magenta", lwd = 3)
lines(x_imc, pred_imc_09, col = "brown", lwd = 2)

legend(
  "topleft",
  legend = c("span = 0.4", "span = 0.75", "span = 0.9"),
  col = c("navy", "magenta", "brown"),
  lwd = c(2, 3, 2)
)

# *****************************************************
# c. MODELO 3: PRESIÓN SISTÓLICA Y RIESGO ALTO CON LOESS
# *****************************************************

loess_presion_04 <- loess(riesgo_num ~ presion_sistolica, data = datos, span = 0.4)
loess_presion_075 <- loess(riesgo_num ~ presion_sistolica, data = datos, span = 0.75)
loess_presion_09 <- loess(riesgo_num ~ presion_sistolica, data = datos, span = 0.9)

x_presion <- seq(min(datos$presion_sistolica), max(datos$presion_sistolica), length.out = 100)

pred_presion_04 <- predict(loess_presion_04, newdata = data.frame(presion_sistolica = x_presion))
pred_presion_075 <- predict(loess_presion_075, newdata = data.frame(presion_sistolica = x_presion))
pred_presion_09 <- predict(loess_presion_09, newdata = data.frame(presion_sistolica = x_presion))

plot(
  datos$presion_sistolica,
  datos$riesgo_num,
  main = "Comparación de span en LOESS: Presión sistólica",
  xlab = "Presión sistólica",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(x_presion, pred_presion_04, col = "navy", lwd = 2)

lines(x_presion, pred_presion_075, col = "yellow", lwd = 3)
lines(x_presion, pred_presion_09, col = "brown", lwd = 2)

legend(
  "topleft",
  legend = c("span = 0.4", "span = 0.75", "span = 0.9"),
  col = c("navy", "yellow", "brown"),
  lwd = c(2, 3, 2)
)

# *****************************************************
# d. MODELO 4: GLUCOSA EN AYUNO Y RIESGO ALTO CON LOESS
# *****************************************************

loess_glucosa_04 <- loess(riesgo_num ~ glucosa_ayuno, data = datos, span = 0.4)
loess_glucosa_075 <- loess(riesgo_num ~ glucosa_ayuno, data = datos, span = 0.75)
loess_glucosa_09 <- loess(riesgo_num ~ glucosa_ayuno, data = datos, span = 0.9)

x_glucosa <- seq(min(datos$glucosa_ayuno), max(datos$glucosa_ayuno), length.out = 100)

pred_glucosa_04 <- predict(loess_glucosa_04, newdata = data.frame(glucosa_ayuno = x_glucosa))
pred_glucosa_075 <- predict(loess_glucosa_075, newdata = data.frame(glucosa_ayuno = x_glucosa))
pred_glucosa_09 <- predict(loess_glucosa_09, newdata = data.frame(glucosa_ayuno = x_glucosa))

plot(
  datos$glucosa_ayuno,
  datos$riesgo_num,
  main = "Comparación de span en LOESS: Glucosa en ayuno",
  xlab = "Glucosa en ayuno",
  ylab = "Riesgo alto: 1 = alto, 0 = bajo",
  pch = 19,
  col = "gray",
  ylim = c(-0.1, 1.1)
)

lines(x_glucosa, pred_glucosa_04, col = "navy", lwd = 2)

lines(x_glucosa, pred_glucosa_075, col = "orange", lwd = 3)
lines(x_glucosa, pred_glucosa_09, col = "brown", lwd = 2)

legend(
  "topleft",
  legend = c("span = 0.4", "span = 0.75", "span = 0.9"),
  col = c("navy", "orange", "brown"),
  lwd = c(2, 3, 2)
)

# ============================================================
# X. REVISIÓN DE ERRORES DE LOS MODELOS LOESS
# ============================================================

# Calculamos las predicciones de cada modelo sobre los mismos datos.
# Esto nos permite comparar lo que el modelo predijo contra el valor real.
pred_edad_datos <- predict(modelo_loess_edad)
pred_imc_datos <- predict(modelo_loess_imc)
pred_presion_datos <- predict(modelo_loess_presion)
pred_glucosa_datos <- predict(modelo_loess_glucosa)

# MSE = promedio de los errores al cuadrado.
mse <- function(real, predicho) {
  mean((real - predicho)^2)
}

# RMSE = raíz cuadrada del MSE.
rmse <- function(real, predicho) {
  sqrt(mean((real - predicho)^2))
}

# Tabla para comparar los errores de cada modelo.
tabla_errores_loess <- data.frame(
  variable = c("edad", "imc", "presion_sistolica", "glucosa_ayuno"),
  MSE = c(
    mse(datos$riesgo_num, pred_edad_datos),
    mse(datos$riesgo_num, pred_imc_datos),
    mse(datos$riesgo_num, pred_presion_datos),
    mse(datos$riesgo_num, pred_glucosa_datos)
  ),
  RMSE = c(
    rmse(datos$riesgo_num, pred_edad_datos),
    rmse(datos$riesgo_num, pred_imc_datos),
    rmse(datos$riesgo_num, pred_presion_datos),
    rmse(datos$riesgo_num, pred_glucosa_datos)
  )
)

# Mostramos la tabla de errores.
tabla_errores_loess

# ============================================================
# IX. GRÁFICAS DE RESIDUOS DE LOS MODELOS LOESS
# ============================================================

# Residuos del modelo con edad.
plot(
  datos$edad,
  modelo_loess_edad$residuals,
  main = "Residuos LOESS: Edad",
  xlab = "Edad",
  ylab = "Residuo",
  pch = 19,
  col = "gray"
)
abline(h = 0, col = "cyan", lty = 2, lwd = 2)

# Residuos del modelo con IMC.
plot(
  datos$imc,
  modelo_loess_imc$residuals,
  main = "Residuos LOESS: IMC",
  xlab = "IMC",
  ylab = "Residuo",
  pch = 19,
  col = "gray"
)
abline(h = 0, col = "magenta", lty = 2, lwd = 2)

# Residuos del modelo con presión sistólica.
plot(
  datos$presion_sistolica,
  modelo_loess_presion$residuals,
  main = "Residuos LOESS: Presión sistólica",
  xlab = "Presión sistólica",
  ylab = "Residuo",
  pch = 19,
  col = "gray"
)
abline(h = 0, col = "yellow", lty = 2, lwd = 2)

# Residuos del modelo con glucosa en ayuno.
plot(
  datos$glucosa_ayuno,
  modelo_loess_glucosa$residuals,
  main = "Residuos LOESS: Glucosa en ayuno",
  xlab = "Glucosa en ayuno",
  ylab = "Residuo",
  pch = 19,
  col = "gray"
)
abline(h = 0, col = "orange", lty = 2, lwd = 2)

# ============================================================
# XII. INCISO 2: MODELO GAM CON TODAS LAS VARIABLES INDEPENDIENTES
# ============================================================

library(mgcv)
datos$riesgo_alto <- factor(datos$riesgo_alto)
datos$sexo <- factor(datos$sexo)
datos$actividad_fisica <- factor(datos$actividad_fisica)
datos$fumador <- factor(datos$fumador)

# Solo para confirmar la variable numérica de riesgo.
# riesgo_num = 1 si la persona tiene riesgo alto.
# riesgo_num = 0 si la persona tiene riesgo bajo.
datos$riesgo_num <- 0
datos$riesgo_num[datos$riesgo_alto == "alto"] <- 1

# Ajusto el modelo GAM.
# Las variables numéricas entran con s(), porque quiero permitir relaciones no lineales suaves.
# Las variables categóricas entran directo como factores.
modelo_gam <- gam(
  riesgo_num ~ 
    s(edad, k = 10) +
    s(imc, k = 10) +
    s(presion_sistolica, k = 10) +
    s(glucosa_ayuno, k = 10) +
    sexo +
    actividad_fisica +
    fumador,
  data = datos,
  family = binomial(), #en lugar de family = gaussian()
  method = "REML" 
)

summary(modelo_gam)

Family: binomial 
Link function: logit 

Formula:
riesgo_num ~ s(edad, k = 10) + s(imc, k = 10) + s(presion_sistolica, 
    k = 10) + s(glucosa_ayuno, k = 10) + sexo + actividad_fisica + 
    fumador

Parametric coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 -8.1698     1.2247  -6.671 2.55e-11 ***
sexoM                        1.0234     0.4998   2.048 0.040593 *  
actividad_fisicamoderado     1.5616     0.7021   2.224 0.026135 *  
actividad_fisicasedentario   2.6142     0.7609   3.436 0.000591 ***
fumadorsi                    3.2808     0.7169   4.576 4.73e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                       edf Ref.df Chi.sq  p-value    
s(edad)              1.000  1.000 54.614  < 2e-16 ***
s(imc)               3.231  4.042 36.038 4.37e-07 ***
s(presion_sistolica) 3.682  4.592 15.658  0.00573 ** 
s(glucosa_ayuno)     1.000  1.000  2.238  0.13472    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.832   Deviance explained = 81.8%
-REML = 62.379  Scale est. = 1         n = 500

# ============================================================
# XIII. GRÁFICAS DE FUNCIONES SUAVES DEL GAM
# ============================================================

# Para ver la contribución parcial de cada variable numérica.
par(mfrow = c(2, 2))

plot(
  modelo_gam,
  shade = TRUE,
  main = "Funciones suaves del modelo GAM"
)

par(mfrow = c(1, 1)) 

# para ver si la flexibilidad elegida para las funciones suaves parece adecuada.
gam.check(modelo_gam)

Method: REML   Optimizer: outer newton
full convergence after 12 iterations.
Gradient range [-1.097051e-05,1.838552e-06]
(score 62.37901 & scale 1).
Hessian positive definite, eigenvalue range [5.504187e-06,0.6052596].
Model rank =  41 / 41 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                       k'  edf k-index p-value
s(edad)              9.00 1.00    1.00    0.41
s(imc)               9.00 3.23    1.01    0.54
s(presion_sistolica) 9.00 3.68    0.95    0.10
s(glucosa_ayuno)     9.00 1.00    0.97    0.23

# ============================================================
# XIV. PREDICCIÓN Y EVALUACIÓN DEL GAM
# ============================================================

# Para mí: type = "response" hace que las predicciones salgan entre 0 y 1.
# Valores cercanos a 1 indican mayor tendencia de riesgo alto.
pred_gam <- predict(
  modelo_gam,
  type = "response"
)

summary(pred_gam)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 8.420e-06 7.048e-03 2.980e-01 7.717e-01 1.000e+00 

# Convertimos las predicciones numéricas en categorías.
# Si la predicción es mayor o igual a 0.5, decimos "alto".
# Si es menor a 0.5, decimos "bajo".
pred_gam_categoria <- ifelse(
  pred_gam >= 0.5,
  "alto",
  "bajo"
)

# Convertimos la predicción a factor.
pred_gam_categoria <- factor(
  pred_gam_categoria,
  levels = c("bajo", "alto")
)

# Nos aseguramos de que la variable real tenga el mismo orden.
datos$riesgo_alto <- factor(
  datos$riesgo_alto,
  levels = c("bajo", "alto")
)

# Hacemos una tabla para comparar la categoría real
# contra la categoría predicha por el modelo.
tabla_prediccion_gam <- table(
  Real = datos$riesgo_alto,
  Predicho = pred_gam_categoria
)

tabla_prediccion_gam

      Predicho
Real   bajo alto
  bajo  343    8
  alto   12  137

# Calculamos el porcentaje total de aciertos.
porcentaje_aciertos_gam <- sum(diag(tabla_prediccion_gam)) / sum(tabla_prediccion_gam)

porcentaje_aciertos_gam

[1] 0.96

# Revisamos porcentajes por fila.
# Esto nos ayuda a ver qué porcentaje de "bajo" y qué porcentaje de "alto"
# fueron clasificados correctamente.
prop.table(tabla_prediccion_gam, 1)

      Predicho
Real         bajo       alto
  bajo 0.97720798 0.02279202
  alto 0.08053691 0.91946309

# ============================================================
# XVI. INCISO 3: ÁRBOL DE DECISIÓN
# ============================================================

library(rpart)

# NOTAS PARA MÍ: 
# 1. method = "class" se usa porque la variable dependiente es categórica.
# 2. maxdepth controla qué tan profundo puede crecer el árbol.
# 3. minsplit indica el mínimo de observaciones necesarias para intentar dividir un nodo.
# 4. cp controla la complejidad del árbol y ayuda a evitar que crezca demasiado.
arbol_decision <- rpart(
  riesgo_alto ~ edad + imc + presion_sistolica + glucosa_ayuno +
    sexo + actividad_fisica + fumador,
  data = datos,
  method = "class",
  control = rpart.control(
    maxdepth = 4,
    minsplit = 10,
    cp = 0.01
  )
)

# Mostramos el árbol en formato de texto.
arbol_decision

n= 500 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 500 149 bajo (0.70200000 0.29800000)  
   2) edad< 65.5 385  43 bajo (0.88831169 0.11168831)  
     4) edad< 51.5 286   5 bajo (0.98251748 0.01748252) *
     5) edad>=51.5 99  38 bajo (0.61616162 0.38383838)  
      10) imc< 28.95 67  15 bajo (0.77611940 0.22388060)  
        20) fumador=no 50   5 bajo (0.90000000 0.10000000) *
        21) fumador=si 17   7 alto (0.41176471 0.58823529) *
      11) imc>=28.95 32   9 alto (0.28125000 0.71875000) *
   3) edad>=65.5 115   9 alto (0.07826087 0.92173913)  
     6) imc< 21.05 11   5 alto (0.45454545 0.54545455)  
      12) actividad_fisica=activo 3   0 bajo (1.00000000 0.00000000) *
      13) actividad_fisica=moderado,sedentario 8   2 alto (0.25000000 0.75000000) *
     7) imc>=21.05 104   4 alto (0.03846154 0.96153846) *

plot(
  arbol_decision,
  uniform = TRUE,
  main = "Árbol de decisión para riesgo cardiovascular",
  margin = 0.1
)

text(
  arbol_decision,
  use.n = TRUE,
  cex = 0.7,
  pretty = 0
)

# ============================================================
# XVII. PREDICCIÓN Y EVALUACIÓN DEL ÁRBOL
# ============================================================

pred_arbol <- predict(
  arbol_decision,
  type = "class"
)

# Creamos una tabla para comparar: 
# filas son las categorías reales
# columnas las categorías predichas.
tabla_arbol <- table(
  Real = datos$riesgo_alto,
  Predicho = pred_arbol
)

tabla_arbol

      Predicho
Real   bajo alto
  bajo  329   22
  alto   10  139

# Para ver el porcentaje total de aciertos.
porcentaje_aciertos_arbol <- sum(diag(tabla_arbol)) / sum(tabla_arbol)

porcentaje_aciertos_arbol

[1] 0.936

# Para ver los aciertos individuales en porcentajes. 
prop.table(tabla_arbol, 1)

      Predicho
Real         bajo       alto
  bajo 0.93732194 0.06267806
  alto 0.06711409 0.93288591

# Para ver qué variables fueron más importantes para construir las divisiones del árbol.
arbol_decision$variable.importance

             edad               imc           fumador presion_sistolica 
      136.2827130        14.8698338         6.0482880         2.8833992 
 actividad_fisica     glucosa_ayuno 
        2.4545455         0.9944251