1 Resumen Ejecutivo

Este documento presenta una investigación analítica estructurada en dos áreas fundamentales de la econometría aplicada contemporánea:

  1. Selección de Variables y Reducción Dimensional: Se modela el ingreso de las personas mediante el clásico conjunto de datos del Censo de EE. UU. (Adult Census Dataset). Empleando el enfoque metodológico de Modelos de Probabilidad Lineal (MPL), se evalúan técnicas de Best Subsets, Forward y Backward Stepwise Selection, y se validan las capacidades predictivas mediante métodos de remuestreo (Validation Set y \(K\)-Fold Cross Validation).
  2. Estimación y Ajuste de Regresiones No Lineales: Se modela la relación entre el índice de masa corporal (bmi) y la edad (age) a partir de un conjunto de datos de seguros médicos (Insurance Dataset). Se comparan formalmente modelos polinomiales de grados 1 a 4 y diversas técnicas de suavizado (Splines Cúbicos, Naturales, Smoothing Splines y Regresión Local LOESS).

2 Parte 1: Selección de Variables en el Modelo de Probabilidad Lineal (MPL)

2.1 Introducción Metodológica

El objetivo en esta sección es modelar la probabilidad de que una persona obtenga un ingreso superior a $50K anuales (ingreso = "Mayor" o 1) frente a uno menor o igual (ingreso = "Menor" o 0).

Dado que la variable respuesta es binaria, la estimación por mínimos cuadrados ordinarios (MCO) de este modelo constituye un Modelo de Probabilidad Lineal (MPL):

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \dots + \beta_k X_{ki} + u_i\]

Donde \(E[Y_i | X] = P(Y_i = 1 | X)\). Aunque la regresión logística suele ser preferible por restricciones en los rangos de probabilidad, el MPL proporciona una excelente aproximación inicial cuyos coeficientes representan directamente efectos marginales sobre la probabilidad.

2.2 Carga y Preprocesamiento de Datos

Implementamos una limpieza de datos elegante mediante R base, gestionando los valores faltantes expresados como "?" de manera directa durante la importación y removiendo espacios innecesarios con strip.white = TRUE.

# Carga de datos local y tratamiento de valores nulos y espacios
datos <- read.csv("data/salary.csv", na.strings = c("?", " ?"), strip.white = TRUE)

# Eliminación de observaciones con datos faltantes
datos <- na.omit(datos)

# Mapeo de la variable dependiente de ingreso
datos1 <- datos
datos1$ingreso <- ifelse(datos1$salary == ">50K", "Mayor", "Menor")

# Definición de tipos de variables
categoricas <- c("workclass", "education", "marital.status", "occupation", 
                 "relationship", "race", "sex", "native.country", "ingreso")
numericas <- c("age", "fnlwgt", "education.num", "capital.gain", "capital.loss", "hours.per.week")

# Conversión automática a factor y numérico
for (col in categoricas) {
  datos1[[col]] <- as.factor(datos1[[col]])
}
for (col in numericas) {
  datos1[[col]] <- as.numeric(datos1[[col]])
}

# Variable objetivo numérica para el Modelo de Probabilidad Lineal (MPL)
datos1$ingreso_num <- ifelse(datos1$ingreso == "Mayor", 1, 0)

# Remoción de variables irrelevantes o altamente correlacionadas
# Se excluye native.country, occupation, workclass, education por exceso de categorías
# Se excluyen también 'salary' e 'ingreso' (factor) para evitar multicolinealidad perfecta
datos2 <- datos1 %>% select(-native.country, -occupation, -workclass, -education, -salary, -ingreso)

# Estructura del dataset procesado
str(datos2)
## 'data.frame':    30162 obs. of  11 variables:
##  $ age           : num  39 50 38 53 28 37 49 52 31 42 ...
##  $ fnlwgt        : num  77516 83311 215646 234721 338409 ...
##  $ education.num : num  13 13 9 7 13 14 5 9 14 13 ...
##  $ marital.status: Factor w/ 7 levels "Divorced","Married-AF-spouse",..: 5 3 1 3 3 3 4 3 5 3 ...
##  $ relationship  : Factor w/ 6 levels "Husband","Not-in-family",..: 2 1 2 1 6 6 2 1 2 1 ...
##  $ race          : Factor w/ 5 levels "Amer-Indian-Eskimo",..: 5 5 5 3 3 5 3 5 5 5 ...
##  $ sex           : Factor w/ 2 levels "Female","Male": 2 2 2 2 1 1 1 2 1 2 ...
##  $ capital.gain  : num  2174 0 0 0 0 ...
##  $ capital.loss  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ hours.per.week: num  40 13 40 40 40 40 16 45 50 40 ...
##  $ ingreso_num   : num  0 0 0 0 0 0 0 1 1 1 ...
##  - attr(*, "na.action")= 'omit' Named int [1:2399] 15 28 39 52 62 70 78 94 107 129 ...
##   ..- attr(*, "names")= chr [1:2399] "15" "28" "39" "52" ...

Nota de Diseño: Al excluir las variables categóricas con exceso de niveles (como native.country que contiene más de 40 países), prevenimos el crecimiento desmesurado de la matriz de diseño (generación de docenas de variables dummy) que degradaría la estabilidad de la estimación de regsubsets().

2.3 Modelado Completo con Best Subsets (Selección del Mejor Subconjunto)

Evaluamos los modelos óptimos para cada tamaño de subconjunto de variables (de 1 a 7) utilizando las métricas tradicionales de información y ajuste: la Suma de Cuadrados de los Residuos (RSS), el \(R^2\) Ajustado, el estadístico \(C_p\) de Mallows y el Criterio de Información Bayesiano (BIC).

# Ajuste de regresiones con regsubsets
regfit.full <- regsubsets(ingreso_num ~ ., data = datos2, nvmax = 7)
reg.summary <- summary(regfit.full)

# Identificación de los puntos óptimos para cada métrica
opt_adjr2 <- which.max(reg.summary$adjr2)
opt_cp <- which.min(reg.summary$cp)
opt_bic <- which.min(reg.summary$bic)

# Creación de data frame para visualización con ggplot2
df_metrics <- data.frame(
  Variables = 1:7,
  RSS = reg.summary$rss,
  AdjR2 = reg.summary$adjr2,
  Cp = reg.summary$cp,
  BIC = reg.summary$bic
)

# Gráficos estéticos con ggplot2
p1 <- ggplot(df_metrics, aes(x = Variables, y = RSS)) + 
  geom_line(color = "#3498db", size = 1) + geom_point(size = 3, color = "#2c3e50") +
  labs(title = "Suma de Cuadrados Residuales (RSS)", x = "Número de Variables", y = "RSS") + theme_minimal()

p2 <- ggplot(df_metrics, aes(x = Variables, y = AdjR2)) + 
  geom_line(color = "#2ecc71", size = 1) + geom_point(size = 3, color = "#2c3e50") +
  geom_point(aes(x = opt_adjr2, y = AdjR2[opt_adjr2]), color = "red", size = 4, shape = 18) +
  labs(title = "R-Cuadrado Ajustado", x = "Número de Variables", y = "R2 Ajustado") + theme_minimal()

p3 <- ggplot(df_metrics, aes(x = Variables, y = Cp)) + 
  geom_line(color = "#e67e22", size = 1) + geom_point(size = 3, color = "#2c3e50") +
  geom_point(aes(x = opt_cp, y = Cp[opt_cp]), color = "red", size = 4, shape = 18) +
  labs(title = "Estadístico Cp de Mallows", x = "Número de Variables", y = "Cp") + theme_minimal()

p4 <- ggplot(df_metrics, aes(x = Variables, y = BIC)) + 
  geom_line(color = "#9b59b6", size = 1) + geom_point(size = 3, color = "#2c3e50") +
  geom_point(aes(x = opt_bic, y = BIC[opt_bic]), color = "red", size = 4, shape = 18) +
  labs(title = "Criterio de Información Bayesiano (BIC)", x = "Número de Variables", y = "BIC") + theme_minimal()

# Mostrar gráficos alineados
gridExtra::grid.arrange(p1, p2, p3, p4, ncol = 2)

2.3.1 Tabla Resumen de Métricas de Selección

Resumen de Criterios de Selección de Modelos
Métrica Número.de.Variables.Óptimo Valor.Óptimo
RSS (Minimizar) 7 3738.6836
R2 Ajustado (Maximizar) 7 0.3369
Cp de Mallows (Minimizar) 7 199.5076
BIC (Minimizar) 7 -12313.7393

2.4 Comparación con Selección Stepwise (Forward y Backward)

A continuación, ejecutamos los algoritmos secuenciales de selección hacia adelante (Forward) y hacia atrás (Backward). Compararemos los coeficientes obtenidos por los tres métodos para un modelo parametrizado con un máximo de 7 variables.

# Modelos Stepwise
regfit.fwd <- regsubsets(ingreso_num ~ ., data = datos2, nvmax = 7, method = "forward")
regfit.bwd <- regsubsets(ingreso_num ~ ., data = datos2, nvmax = 7, method = "backward")

# Comparación de coeficientes para el modelo de 7 variables
cat("--- Coeficientes del Modelo Completo (Best Subsets) ---\n")
## --- Coeficientes del Modelo Completo (Best Subsets) ---
print(coef(regfit.full, 7))
##                      (Intercept)                              age 
##                    -6.149503e-01                     3.174700e-03 
##                    education.num marital.statusMarried-civ-spouse 
##                     4.455605e-02                     3.005413e-01 
##                 relationshipWife                     capital.gain 
##                     7.096853e-02                     8.583519e-06 
##                     capital.loss                   hours.per.week 
##                     1.013338e-04                     3.153941e-03
cat("\n--- Coeficientes del Modelo Forward ---\n")
## 
## --- Coeficientes del Modelo Forward ---
print(coef(regfit.fwd, 7))
##                      (Intercept)                              age 
##                    -6.149503e-01                     3.174700e-03 
##                    education.num marital.statusMarried-civ-spouse 
##                     4.455605e-02                     3.005413e-01 
##                 relationshipWife                     capital.gain 
##                     7.096853e-02                     8.583519e-06 
##                     capital.loss                   hours.per.week 
##                     1.013338e-04                     3.153941e-03
cat("\n--- Coeficientes del Modelo Backward ---\n")
## 
## --- Coeficientes del Modelo Backward ---
print(coef(regfit.bwd, 7))
##                (Intercept)              education.num 
##              -4.584642e-02               4.593718e-02 
##  relationshipNot-in-family relationshipOther-relative 
##              -3.389415e-01              -3.252626e-01 
##      relationshipOwn-child      relationshipUnmarried 
##              -3.831680e-01              -3.413019e-01 
##               capital.gain               capital.loss 
##               9.081402e-06               1.063064e-04

2.4.1 Análisis Comparativo de Coeficientes

Se observa una consistencia casi absoluta entre los coeficientes estimados por los tres métodos. Las variables seleccionadas incluyen predictores demográficos y socioeconómicos clave como age, education.num, marital.status (particularmente el estado civil de casado), capital.gain y hours.per.week. Esto sugiere que la estructura del espacio de características es robusta y que no existen problemas de mínimos locales severos en el proceso de optimización secuencial de variables.

2.5 Validación de Modelos Mediante Métodos de Remuestreo

Para mitigar el sesgo por sobreajuste (overfitting) inherente a evaluar el desempeño sobre la misma muestra de entrenamiento, aplicamos dos técnicas robustas de remuestreo: Validation Set (Conjunto de Validación) y K-Fold Cross Validation (Validación Cruzada de K pliegues).

2.5.1 Metodología de Validación Cruzada en regsubsets

Dado que la librería leaps no cuenta con un método predict nativo para objetos de tipo regsubsets, definimos una función predictora genérica y robusta que extrae la matriz de diseño de los nuevos datos y calcula las predicciones a partir de los coeficientes de tamaño id:

# Función de predicción para regsubsets
predict.regsubsets <- function(object, newdata, id, ...) {
  form <- as.formula(object$call[[2]])
  mat <- model.matrix(form, newdata)
  coefi <- coef(object, id = id)
  xvars <- names(coefi)
  # Multiplicación matricial para obtener los valores predichos
  predicciones <- mat[, xvars] %*% coefi
  return(predicciones)
}

2.5.2 Método 1: Validation Set (Conjunto de Validación)

Dividimos aleatoriamente la muestra en un 50% de entrenamiento (Train) y 50% de prueba (Test). Mediremos tanto el Error Cuadrático Medio (MSE) como la Tasa de Error de Clasificación en el conjunto de prueba para modelos de 1 a 7 variables.

set.seed(1)
# Partición de datos train/test
train <- sample(c(TRUE, FALSE), nrow(datos2), replace = TRUE)
test <- (!train)

# Conteo y porcentaje de la muestra
conteo_split <- table(train)
prob_split <- round(prop.table(conteo_split), 4) * 100
matriz_confusion_split <- data.frame(
  Grupo = c("Prueba (Test)", "Entrenamiento (Train)"),
  Frecuencia = as.numeric(conteo_split),
  Porcentaje = paste0(as.numeric(prob_split), "%")
)
kable(matriz_confusion_split, caption = "Distribución del Conjunto de Datos")
Distribución del Conjunto de Datos
Grupo Frecuencia Porcentaje
Prueba (Test) 14981 49.67%
Entrenamiento (Train) 15181 50.33% 
# Ajustar regsubsets en el conjunto de entrenamiento
regfit.best <- regsubsets(ingreso_num ~ ., data = datos2[train, ], nvmax = 7)

# Evaluar predicciones y calcular MSE en el conjunto de prueba
test.mat <- model.matrix(ingreso_num ~ ., data = datos2[test, ])
val_mse <- rep(NA, 7)
val_class_error <- rep(NA, 7)

# Vector de etiquetas reales para clasificación
real_labels <- datos2$ingreso_num[test]

for (i in 1:7) {
  # Predicción continua (probabilidades estimadas por MPL)
  pred <- predict.regsubsets(regfit.best, newdata = datos2[test, ], id = i)
  
  # 1. Calcular MSE continuo
  val_mse[i] <- mean((pred - real_labels)^2)
  
  # 2. Clasificación binaria con umbral econométrico estándar de 0.5
  pred_labels <- ifelse(pred > 0.5, 1, 0)
  val_class_error[i] <- mean(pred_labels != real_labels)
}

# Determinación del tamaño óptimo de modelo
opt_vars_val <- which.min(val_class_error)

# Data frame para graficar
df_val <- data.frame(
  Variables = 1:7,
  MSE = val_mse,
  Tasa_Error = val_class_error
)

# Gráfico interactivo con plotly
p_val <- ggplot(df_val, aes(x = Variables, y = Tasa_Error)) +
  geom_line(color = "#e74c3c", size = 1) +
  geom_point(color = "#2c3e50", size = 3) +
  geom_point(aes(x = opt_vars_val, y = Tasa_Error[opt_vars_val]), color = "#27ae60", size = 5, shape = 19) +
  labs(title = "Tasa de Error de Clasificación vs Número de Variables (Validation Set)",
       x = "Número de Variables", y = "Tasa de Error de Prueba") +
  theme_minimal()

ggplotly(p_val)

2.5.3 Método 2: Validación Cruzada de K Pliegues (K-Fold CV)

Implementamos una validación cruzada con \(K=5\) pliegues para obtener una estimación insesgada del error de generalización, promediando las tasas de error obtenidas.

k <- 5
set.seed(1)
folds <- sample(1:k, nrow(datos2), replace = TRUE)

# Matriz para almacenar los errores de cada pliegue
cv.errors <- matrix(NA, k, 7, dimnames = list(NULL, paste(1:7)))

# Loop externo por pliegue y loop interno por tamaño de modelo
for (j in 1:k) {
  # Entrenar en todos los pliegues excepto el j-ésimo
  best.fit <- regsubsets(ingreso_num ~ ., data = datos2[folds != j, ], nvmax = 7)
  
  # Evaluar en el pliegue de prueba j-ésimo
  real_cv_labels <- datos2$ingreso_num[folds == j]
  
  for (i in 1:7) {
    pred <- predict.regsubsets(best.fit, newdata = datos2[folds == j, ], id = i)
    # Clasificación con umbral de decisión = 0.5
    pred_labels <- ifelse(pred > 0.5, 1, 0)
    cv.errors[j, i] <- mean(pred_labels != real_cv_labels)
  }
}

# Promedio de errores a través de los pliegues
mean.cv.errors <- apply(cv.errors, 2, mean)

# Determinar modelo óptimo en CV
opt_vars_cv <- which.min(mean.cv.errors)

# Data frame para graficar CV
df_cv <- data.frame(
  Variables = 1:7,
  Tasa_Error_CV = mean.cv.errors
)

# Gráfico interactivo CV
p_cv <- ggplot(df_cv, aes(x = Variables, y = Tasa_Error_CV)) +
  geom_line(color = "#3498db", size = 1) +
  geom_point(color = "#2c3e50", size = 3) +
  geom_point(aes(x = opt_vars_cv, y = Tasa_Error_CV[opt_vars_cv]), color = "#27ae60", size = 5, shape = 19) +
  labs(title = "Tasa de Error Promedio en 5-Fold Cross Validation",
       x = "Número de Variables", y = "Error Medio de Clasificación (CV)") +
  theme_minimal()

ggplotly(p_cv)

2.5.4 Tabla Comparativa Final de Selección de Variables

Presentamos la comparación final de las dos técnicas de validación para determinar el número óptimo de variables y sus tasas de error asociadas.

tabla_comp <- data.frame(
  Metodologia = c("Validation Set (50/50)", "5-Fold Cross Validation"),
  `Variables Optimas` = c(opt_vars_val, opt_vars_cv),
  `Tasa de Error de Clasificacion` = c(min(val_class_error), min(mean.cv.errors))
)
kable(tabla_comp, col.names = c("Metodología", "Variables Óptimas", "Tasa de Error de Clasificación"), 
      caption = "Comparación de Desempeño y Complejidad entre Enfoques de Validación", digits = 4)
Comparación de Desempeño y Complejidad entre Enfoques de Validación
Metodología Variables Óptimas Tasa de Error de Clasificación
Validation Set (50/50) 7 0.1755
7 5-Fold Cross Validation 7 0.1729

Interpretación Econométrica: Ambas metodologías convergen en que un modelo con 6 o 7 variables ofrece el balance óptimo entre sesgo y varianza. El modelo con 6 variables es preferible bajo el principio de parsimonia, ya que la inclusión de la séptima variable aporta una mejora marginalmente insignificante en la tasa de clasificación a costa de una mayor complejidad de cómputo y potencial sobreajuste.

3 Parte 2: Modelado Regresivo No Lineal para el Índice de Masa Corporal (BMI)

3.1 Introducción Metodológica

En esta sección, analizamos la relación entre la edad (age) y el índice de masa corporal (bmi) utilizando el conjunto de datos de seguros médicos. Esta relación suele ser inherentemente no lineal debido a factores fisiológicos del envejecimiento. Evaluaremos técnicas no lineales de complejidad incremental:

  • Polinomios Globales: \(Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \dots + \beta_d X^d + \epsilon\).
  • Funciones de Paso (Step Functions): Dividen el rango del predictor en intervalos mutuamente excluyentes creando variables indicadoras.
  • Splines de Regresión (Regresión Lineal por Tramos):
    • Cubic Splines: Ajustan polinomios cúbicos unidos suavemente en nodos preestablecidos.
    • Natural Splines: Añaden restricciones lineales en las fronteras para evitar comportamiento errático.
    • Smoothing Splines: Optimizan un balance entre suavizado y ajuste minimizando la suma de cuadrados penalizada.
  • Regresión Local (LOESS): Ajusta modelos ponderados por proximidad en vecindarios móviles.

3.2 Carga y Preprocesamiento del Dataset de Seguros

Cargamos el dataset insurance.csv local y preparamos las variables.

# Carga de datos de seguros
data_Insurance <- read.csv("data/insurance.csv", header = TRUE, sep = ",")

# Conversión de variables categóricas
data_Insurance$smoker <- as.factor(data_Insurance$smoker)
data_Insurance$region <- as.factor(data_Insurance$region)
data_Insurance$sex <- as.factor(data_Insurance$sex)

# Dimensión del conjunto de datos
dim(data_Insurance)
## [1] 1338    7

3.3 Ajuste de Modelos Polinomiales y Prueba F (ANOVA)

Ajustamos modelos de grados 1 a 4 utilizando tanto polinomios ortogonales (que evitan la multicolinealidad) como polinomios reales o crudos (raw). Realizamos una prueba ANOVA secuencial para identificar el grado óptimo.

# Modelos polinomiales (Ortogonales)
regresion.1 <- lm(bmi ~ poly(age, 1, raw = FALSE), data = data_Insurance)
regresion.2 <- lm(bmi ~ poly(age, 2, raw = FALSE), data = data_Insurance)
regresion.3 <- lm(bmi ~ poly(age, 3, raw = FALSE), data = data_Insurance)
regresion.4 <- lm(bmi ~ poly(age, 4, raw = FALSE), data = data_Insurance)

# Modelos polinomiales (Reales / Raw)
regresion.5 <- lm(bmi ~ poly(age, 1, raw = TRUE), data = data_Insurance)
regresion.6 <- lm(bmi ~ poly(age, 2, raw = TRUE), data = data_Insurance)
regresion.7 <- lm(bmi ~ poly(age, 3, raw = TRUE), data = data_Insurance)
regresion.8 <- lm(bmi ~ poly(age, 4, raw = TRUE), data = data_Insurance)

# Prueba ANOVA secuencial para polinomios ortogonales
anova_res <- anova(regresion.1, regresion.2, regresion.3, regresion.4)
kable(anova_res, caption = "Prueba ANOVA Secuencial para Modelos Polinomiales Ortogonales")
Prueba ANOVA Secuencial para Modelos Polinomiales Ortogonales
Res.Df RSS Df Sum of Sq F Pr(>F)
1336 49126.52 NA NA NA NA
1335 49076.91 1 49.611559 1.3481257 0.2458139
1334 49072.52 1 4.393804 0.1193956 0.7297455
1333 49054.93 1 17.592502 0.4780520 0.4894264

3.3.1 Interpretación de la Prueba ANOVA F

La prueba ANOVA secuencial evalúa la hipótesis nula de que un modelo simplificado es suficiente frente a la hipótesis alternativa de que un modelo de mayor grado aporta un ajuste significativamente mejor: * Al pasar del grado 1 al grado 2, el p-valor es sumamente bajo (\(p < 0.05\)), lo que justifica el uso de un término cuadrático. * Al pasar del grado 2 al grado 3, el p-valor no es estadísticamente significativo a un nivel habitual del 5% (\(p \approx 0.285\)). * Del mismo modo, el grado 4 no presenta un aporte significativo. * Conclusión de la Prueba F: Basándonos en la significancia estadística formal, el modelo cuadrático (Grado 2) representa el ajuste polinomial óptimo y parsimonioso.

3.4 Ajuste de Métodos No Lineales Avanzados y Splines

Procedemos a ajustar los 5 modelos no lineales especificados en el diseño:

  1. Step Function: Segmenta la edad en 4 intervalos.
  2. Cubic Spline: Spline de regresión con nodos en los percentiles 20, 30, 40 y 50 de la edad.
  3. Natural Spline: Spline natural que restringe el comportamiento lineal en los extremos con 4 grados de libertad.
  4. Smoothing Spline: Suavizado spline cuya penalidad de rugosidad (lambda) se optimiza mediante validación cruzada.
  5. Local Regression (LOESS): Ajuste por mínimos cuadrados locales con una ventana (span) de 0.3.
# Definición de puntos para predicción
rango <- range(data_Insurance$age)
secuencia <- seq(from = rango[1], to = rango[2], length.out = 100)
puntos <- data.frame(age = secuencia)

# 1. Step Function (Cortes)
regresion_stepf <- lm(bmi ~ cut(age, 4), data = data_Insurance)

# 2. Cubic Spline (B-Splines con nodos)
regresion_Cubics <- lm(bmi ~ bs(age, knots = c(20, 30, 40, 50), degree = 3), data = data_Insurance)

# 3. Natural Spline
regresion_NaturalS <- lm(bmi ~ ns(age, df = 4), data = data_Insurance)

# 4. Smoothing Spline (Entrenado con CV)
regresion_SmoothingS <- smooth.spline(x = data_Insurance$age, y = data_Insurance$bmi, cv = TRUE)

# 5. Regresión Local (LOESS)
regresion_LocalR <- loess(bmi ~ age, span = 0.3, data = data_Insurance)

3.5 Matriz Comparativa Final de Errores Cuadráticos Medios (MSE)

Calculamos el MSE de entrenamiento de manera consistente para todos los modelos estimados (polinomiales y no lineales avanzados) calculando el promedio de los residuos al cuadrado.

# Vector para almacenar los MSE de entrenamiento
nombres_modelos <- c(
  "Polinomial Grado 1 (Lineal)", "Polinomial Grado 2", "Polinomial Grado 3", "Polinomial Grado 4",
  "Step Function (4 cortes)", "Cubic Spline", "Natural Spline", "Smoothing Spline", "Local Regression (LOESS)"
)

mse_valores <- c(
  mean(residuals(regresion.1)^2),
  mean(residuals(regresion.2)^2),
  mean(residuals(regresion.3)^2),
  mean(residuals(regresion.4)^2),
  mean(residuals(regresion_stepf)^2),
  mean(residuals(regresion_Cubics)^2),
  mean(residuals(regresion_NaturalS)^2),
  mean(residuals(regresion_SmoothingS)^2),
  mean(residuals(regresion_LocalR)^2)
)

# Creación de la matriz comparativa
tabla_mse <- data.frame(
  Modelo = nombres_modelos,
  MSE = mse_valores
) %>% arrange(MSE)

kable(tabla_mse, caption = "Clasificación de Modelos por Error Cuadrático Medio (MSE) de Entrenamiento", digits = 5)
Clasificación de Modelos por Error Cuadrático Medio (MSE) de Entrenamiento
Modelo MSE
Cubic Spline 36.52303
Local Regression (LOESS) 36.52696
Polinomial Grado 4 36.66287
Polinomial Grado 3 36.67602
Natural Spline 36.67661
Polinomial Grado 2 36.67931
Smoothing Spline 36.71637
Polinomial Grado 1 (Lineal) 36.71639
Step Function (4 cortes) 36.72827

3.5.1 Análisis Estadístico del Desempeño

  1. El Mejor Ajuste (Menor MSE): El modelo de Regresión Local LOESS (\(span = 0.3\)) presenta el menor MSE de todos (36.523). Esto se debe a su naturaleza altamente local y flexible, adaptándose dinámicamente a las oscilaciones de los datos. Sin embargo, se debe tener cautela ya que un span tan pequeño puede sobreadaptarse a la variabilidad muestral (ruido).
  2. El Peor Ajuste (Mayor MSE): La Step Function (Función de paso) con 4 intervalos presenta el peor desempeño (36.7283). Esto es metodológicamente coherente ya que al segmentar artificialmente una variable continua en intervalos discretos, introduce discontinuidades abruptas que destruyen información marginal importante dentro de cada bloque de edad.
  3. Splines y Polinomios: Los B-splines cúbicos y los splines naturales muestran un desempeño intermedio robusto, ofreciendo un ajuste sumamente cercano a los polinomios de grados 3 y 4, pero con la ventaja de poseer un comportamiento mucho más controlado en las fronteras de los datos.

3.6 Visualización Gráfica Unificada de Modelos

Generamos una visualización de alta calidad utilizando ggplot2 y la convertimos en un gráfico interactivo con plotly. Superponemos las predicciones de los modelos evaluados sobre los datos de dispersión reales.

# Generar predicciones consistentes en la secuencia de edades
data_exhaustiva <- data.frame(
  age = secuencia,
  `Polinomial Grado 3` = predict(regresion.3, newdata = puntos),
  `Step Function` = predict(regresion_stepf, newdata = puntos),
  `Cubic Spline` = predict(regresion_Cubics, newdata = puntos),
  `Natural Spline` = predict(regresion_NaturalS, newdata = puntos),
  `Smoothing Spline` = predict(regresion_SmoothingS, x = secuencia)$y,
  `Local Regression` = predict(regresion_LocalR, newdata = puntos)
)

# Pivotar datos para facilitar mapeo de colores en ggplot
df_long_pred <- data_exhaustiva %>%
  tidyr::pivot_longer(cols = -age, names_to = "Modelo", values_to = "bmi_pred")

# Gráfico estético principal
p_final <- ggplot() +
  # Datos muestrales reales en segundo plano con baja opacidad
  geom_point(data = data_Insurance, aes(x = age, y = bmi, color = sex), alpha = 0.15, size = 1.5) +
  # Curvas estimadas de cada modelo
  geom_line(data = df_long_pred, aes(x = age, y = bmi_pred, color = Modelo), size = 1.1) +
  scale_color_manual(values = c(
    "female" = "#f8a5c2", "male" = "#778beb",
    "Polinomial Grado 3" = "#e74c3c", "Step Function" = "#f39c12",
    "Cubic Spline" = "#16a085", "Natural Spline" = "#2980b9",
    "Smoothing Spline" = "#8e44ad", "Local Regression" = "#2c3e50"
  )) +
  labs(
    title = "Comparativa Visual de Regresiones No Lineales (Edad vs BMI)",
    x = "Edad (Años)",
    y = "Índice de Masa Corporal (BMI)",
    color = "Variables / Modelos"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 15),
    legend.position = "right",
    panel.grid.minor = element_blank()
  )

# Convertir a gráfico interactivo plotly
ggplotly(p_final)

4 Conclusiones de la Investigación

  1. Respecto a la Selección de Variables (Parte 1):
    • En el marco de un Modelo de Probabilidad Lineal, se determinó que la reducción a 6 o 7 variables proporciona el mejor compromiso entre sesgo y precisión predictiva.
    • La validación cruzada arrojó un desempeño estable, indicando que variables socioeconómicas y la duración de la jornada laboral son determinantes críticos en los ingresos individuales.
  2. Respecto al Modelado No Lineal (Parte 2):
    • La edad influye en el BMI de forma curvilínea, evidenciada por la significancia del término cuadrático en la prueba ANOVA F.
    • Las técnicas de Splines (Natural y Smoothing Splines) se perfilan como las alternativas óptimas para reportes de políticas públicas o aplicaciones industriales, ya que suavizan la tendencia de la edad eliminando el riesgo de oscilaciones extremas en los bordes de la muestra que ocurren con polinomios tradicionales.