1 Módulo 2 — Regresión lineal múltiple con enfoque predictivo

1.1 Contexto biológico

Usaremos datos morfométricos de tortugas pintadas (Chrysemys picta), un modelo clásico en ecología de reptiles. El objetivo es predecir la longitud de la concha a partir de variables morfológicas y ambientales.

La diferencia entre enfoque explicativo y predictivo es central en bioestadística moderna: en el explicativo nos interesan los coeficientes e inferencia; en el predictivo nos interesa el desempeño en datos nuevos no vistos.

1.2 Datos

# ── Datos morfométricos basados en relaciones alométricas reales ──────────────
# Referencia: Janzen et al. (1992) Ecology
set.seed(42)
n <- 280

datos_tortugas <- tibble(
  longitud_concha_mm = rnorm(n, mean = 138, sd = 22),
  sexo               = sample(c("Hembra", "Macho"), n, replace = TRUE,
                              prob = c(0.55, 0.45))
) %>%
  mutate(
    longitud_concha_mm = longitud_concha_mm + ifelse(sexo == "Hembra", 8, 0),
    peso_g          = 0.0018 * longitud_concha_mm^2.7 +
      ifelse(sexo == "Hembra", 60, 0) + rnorm(n, 0, 35),
    ancho_concha_mm = 0.82 * longitud_concha_mm + rnorm(n, 0, 8),
    altura_concha_mm= 0.45 * longitud_concha_mm + rnorm(n, 0, 6),
    largo_cabeza_mm = 0.28 * longitud_concha_mm + rnorm(n, 0, 5),
    temp_agua_c     = runif(n, 12, 28),
    profundidad_m   = runif(n, 0.2, 3.5),
    edad_años       = pmax(round(longitud_concha_mm/18 + rnorm(n, 0, 1.5)), 1)
  )

glimpse(datos_tortugas)
## Rows: 280
## Columns: 9
## $ longitud_concha_mm <dbl> 168.16109, 133.57664, 153.98883, 159.92298, 154.893…
## $ sexo               <chr> "Macho", "Hembra", "Hembra", "Hembra", "Hembra", "H…
## $ peso_g             <dbl> 1867.4922, 997.0991, 1513.8624, 1645.5372, 1564.640…
## $ ancho_concha_mm    <dbl> 140.24963, 112.67478, 118.26409, 128.53102, 118.946…
## $ altura_concha_mm   <dbl> 80.23671, 57.76213, 75.35438, 68.20350, 59.89479, 7…
## $ largo_cabeza_mm    <dbl> 56.97275, 34.32343, 50.71398, 42.61396, 38.21520, 3…
## $ temp_agua_c        <dbl> 13.69241, 24.35050, 22.17509, 21.47685, 12.18190, 1…
## $ profundidad_m      <dbl> 2.9048596, 0.7727506, 0.5567349, 1.5877060, 0.42013…
## $ edad_años          <dbl> 10, 6, 7, 8, 9, 8, 8, 10, 11, 9, 9, 11, 8, 9, 9, 6,…

Este proceso se realizo para la generación de los datos que serian posteriormente usados para el estudio de las tortugas, dando por resultado un total de 280 tortugas divididas entre machos y hembras, cada una con los valores de sus datos ya definidos.

1.3 Análisis exploratorio

p1 <- ggplot(datos_tortugas, aes(x = longitud_concha_mm, fill = sexo)) +
  geom_histogram(bins = 25, color = "white", alpha = 0.85) +
  scale_fill_manual(values = c("#FF9EBC", "#00BCD4")) +
  labs(title = "Variable respuesta: longitud de concha",
       x = "Longitud (mm)", y = "Frecuencia") + theme_minimal()

p2 <- ggplot(datos_tortugas, aes(x = peso_g, y = longitud_concha_mm,
                                 color = sexo)) +
  geom_point(alpha = 0.45, size = 1.6) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 1) +
  scale_color_manual(values = c("#FF9EBC", "#00BCD4")) +
  labs(title = "Longitud vs Peso", x = "Peso (g)",
       y = "Longitud concha (mm)") + theme_minimal()

p3 <- ggplot(datos_tortugas, aes(x = temp_agua_c, y = longitud_concha_mm,
                                 color = sexo)) +
  geom_point(alpha = 0.45, size = 1.6) +
  geom_smooth(method = "loess", se = TRUE, color = "gray30", linewidth = 1) +
  scale_color_manual(values = c("#FF9EBC", "#00BCD4")) +
  labs(title = "Longitud vs Temperatura", x = "Temperatura (°C)",
       y = "Longitud concha (mm)") + theme_minimal()

p4 <- ggplot(datos_tortugas, aes(x = sexo, y = longitud_concha_mm,
                                 fill = sexo)) +
  geom_violin(alpha = 0.35) +
  geom_boxplot(width = 0.18, alpha = 0.85) +
  scale_fill_manual(values = c("#FF9EBC", "#00BCD4"), guide = "none") +
  labs(title = "Longitud por sexo", x = NULL,
       y = "Longitud concha (mm)") + theme_minimal()

grid.arrange(p1, p2, p3, p4, ncol = 2)

Los gráficos exploratorios revelan que las hembras tienden a poseer caparazones más largos que los machos en promedio. Esto es consistente con la biología reproductiva de Chrysemys picta: dado que el caparazón no puede expandirse para dar cabida a los ovarios y embriones en crecimiento, la selección natural favorece un mayor tamaño corporal en las hembras como medio para maximizar el tamaño de la nidada o el tamaño de los huevos (Litzgus & Smith, 2010). La relación entre peso y longitud es fuertemente positiva en ambos sexos, mientras que la temperatura del agua no muestra ninguna tendencia sistemática con el tamaño corporal.

vars_num <- datos_tortugas %>%
  select(longitud_concha_mm, peso_g, ancho_concha_mm, altura_concha_mm,
         largo_cabeza_mm, temp_agua_c, profundidad_m, edad_años)

corrplot(cor(vars_num, use = "complete.obs"),
         method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.72,
         tl.cex = 0.9, tl.col = "gray20",
         col = COL2("RdBu", 10),
         title = "Correlaciones entre variables morfométricas",
         mar = c(0, 0, 2, 0))

La matriz de correlaciones confirma que los aspectos morfológicos poseen una alta correlación entre sí. La variable con mayor relación con la longitud de la concha es el peso, seguido del ancho, la altura, el largo de cabeza y finalmente la edad, que presenta una asociación considerablemente menor. Esto se explica porque el tamaño corporal y las estructuras anatómicas de las tortugas escalan casi isométricamente. Esto a su vez, se centra en las bases evolutivas de este grupo de animales, siendo que la evolución del caparazón limitó las proporciones corporales del clado, observándose desviaciones importantes solo en grupos ecológicos altamente especializados (Hermanson & Evers, 2024). Las variables ambientales temp_agua_c y profundidad_m muestran correlaciones cercanas a cero, siendo claramente irrelevantes al ser comparadas con todas las variables morfométricas, lo que sugiere que no aportan un valor predictivo real al modelo. Debido a esto, sería valioso reemplazarlas por variables con mayor relevancia biológica, como por ejemplo características del hábitat o variables geoespaciales que capturen la variación ambiental entre sitios de nacimiento, y a su vez nos puedan dar pistas de qué otros posibles factores están afectando el desarrollo de estos animales. ## División train / test y preprocesamiento

set.seed(123)

split <- initial_split(datos_tortugas, prop = 0.75, strata = sexo)
train <- training(split);  test <- testing(split)

cat(sprintf("Entrenamiento: %d | Prueba: %d\n", nrow(train), nrow(test)))
## Entrenamiento: 210 | Prueba: 70
receta <- recipe(longitud_concha_mm ~ ., data = train) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_normalize(all_numeric_predictors()) %>%
  step_corr(all_numeric_predictors(), threshold = 0.90)

receta_prep <- prep(receta, training = train)
train_proc  <- bake(receta_prep, new_data = NULL)
test_proc   <- bake(receta_prep, new_data = test)

cat("Predictores tras preproceso:",
    paste(setdiff(names(train_proc), "longitud_concha_mm"), collapse = ", "), "\n")
## Predictores tras preproceso: ancho_concha_mm, altura_concha_mm, largo_cabeza_mm, temp_agua_c, profundidad_m, edad_años, sexo_Macho

Aqui se distributeronn los datos para alimentar el modelo en entrenamiento y prueba, mientras que tambien se evaluaron las variables a tener en cuenta para realizar las predicciones. ## Ajuste y evaluación del modelo

modelo_lm <- lm(longitud_concha_mm ~ ., data = train_proc)

tidy(modelo_lm) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                         p.value < 0.05  ~ "*",   TRUE ~ "")) %>%
  kable(caption = "Coeficientes del modelo (datos de entrenamiento)",
        col.names = c("Término","Estimado","Error est.","t","p-valor","")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Coeficientes del modelo (datos de entrenamiento)
Término Estimado Error est. t p-valor
(Intercept) 141.5412 0.4492 315.1075 0.0000 ***
ancho_concha_mm 11.4196 0.8743 13.0613 0.0000 ***
altura_concha_mm 5.8791 0.7724 7.6118 0.0000 ***
largo_cabeza_mm 4.3439 0.6843 6.3480 0.0000 ***
temp_agua_c -0.2816 0.4520 -0.6230 0.5340
profundidad_m -0.7981 0.4529 -1.7620 0.0796
edad_años 1.8540 0.6167 3.0064 0.0030 **
sexo_Macho -0.8371 0.4593 -1.8228 0.0698

La tabla de coeficientes permite profundizar la información obtenida mediante la interpretación de las gráficas exploratorias, señalando que, aunque las hembras tienen caparazones más largos en promedio, el sexo no muestra una significancia estadística relevante sobre la longitud del caparazón una vez que se controlan los predictores morfométricos, lo que indica que la diferencia entre sexos es una consecuencia de las relaciones isométricas de tamaño corporal más que un efecto anatómico o fisiológico independiente. Al mismo tiempo, esta tabla nos permite evidenciar cuáles son los predictores con mayor influencia dentro de las predicciones del modelo, siendo que los datos morfométricos representan la información más importante a la hora de evaluar la posible longitud de la concha, lo cual es biológicamente esperable dado que todas estas variables están gobernadas por las mismas relaciones alométricas.

glance(modelo_lm) %>%
  select(r.squared, adj.r.squared, sigma, AIC, BIC) %>%
  mutate(across(everything(), ~ round(.x, 4))) %>%
  kable(caption = "Bondad de ajuste — datos de entrenamiento",
        col.names = c("R²","R² ajustado","RMSE (σ)","AIC","BIC")) %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)
Bondad de ajuste — datos de entrenamiento
R² R² ajustado RMSE (σ) AIC BIC
0.9174 0.9145 6.5093 1392.555 1422.678 
pred_vals <- predict(modelo_lm, newdata = test_proc)

res_test <- tibble(
  observado = test$longitud_concha_mm,
  predicho  = pred_vals,
  residuo   = observado - predicho
)

rmse_v <- sqrt(mean(res_test$residuo^2))
mae_v  <- mean(abs(res_test$residuo))
r2_v   <- cor(res_test$observado, res_test$predicho)^2

tibble(Métrica = c("RMSE (mm)","MAE (mm)","R²"),
       `Prueba (no visto)` = c(round(rmse_v,3), round(mae_v,3), round(r2_v,4))) %>%
  kable(caption = "Métricas de predicción en datos de PRUEBA") %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)
Métricas de predicción en datos de PRUEBA
Métrica Prueba (no visto)
RMSE (mm) 6.3460
MAE (mm) 5.1580
R² 0.8935

Seguido de la evaluación de coeficientes y su relevancia para el funcionamiento del modelo, se decidió explorar más a fondo los datos brindados por este tanto en el periodo de entrenamiento como en el de prueba. Siguiendo esta idea, al comparar el R² entre entrenamiento y prueba se observó una disminución moderada entre la calidad del modelo en ambas fases, siendo que el modelo logró explicar aproximadamente un 92% de la variabilidad durante el entrenamiento, mientras que en la prueba este porcentaje se redujo a un 89%, marcando una diferencia de solo 3 puntos porcentuales que no constituye evidencia preocupante de sobreajuste. Aún así, se debe recalcar que la RMSE obtenida de los datos de prueba es más informativa y confiable que la de entrenamiento, ya que esta primera fue calculada sobre valores completamente desconocidos para el modelo, eliminando el posible sesgo optimista que hubiera podido tener el modelo en caso de un sobreajuste.

ggplot(res_test, aes(x = observado, y = predicho)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed",
              color = "gray50", linewidth = 0.9) +
  geom_point(aes(color = abs(residuo)), size = 2, alpha = 0.7) +
  scale_color_viridis_c(option = "plasma", name = "|Residuo|\n(mm)") +
  annotate("text",
           x = min(res_test$observado) + 3,
           y = max(res_test$predicho) - 3,
           label = sprintf("RMSE = %.2f mm\nR² = %.3f", rmse_v, r2_v),
           hjust = 0, size = 4, color = "gray20") +
  labs(title    = "Observado vs predicho — datos de prueba",
       subtitle = "Línea punteada = predicción perfecta",
       x = "Longitud observada (mm)", y = "Longitud predicha (mm)") +
  theme_minimal(base_size = 12)

El gráfico de observado vs predicho respalda visualmente las conclusiones previas, demostrando que los puntos del conjunto de prueba siguen de cerca la línea de predicción perfecta sin recaer en un sobreajuste que hiciera que los datos perdieran sentido predictivo, confirmando que el modelo generaliza de forma coherente a datos no vistos. ## Validación cruzada 10-fold

set.seed(456)

cv_res <- workflow() %>%
  add_formula(longitud_concha_mm ~ .) %>%
  add_model(linear_reg() %>% set_engine("lm")) %>%
  fit_resamples(
    resamples = vfold_cv(train_proc, v = 10),
    metrics   = metric_set(rmse, rsq, mae)
  )

collect_metrics(cv_res) %>%
  select(.metric, mean, std_err, n) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Validación cruzada 10-fold",
        col.names = c("Métrica","Media","Error est.","Folds")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Validación cruzada 10-fold
Métrica Media Error est. Folds
mae 5.2584 0.2673 10
rmse 6.5165 0.3358 10
rsq 0.9116 0.0101 10

Como una evidencia aún más sólida de la viabilidad del modelo, se realizó una validación cruzada, la cual terminó por complementar la evaluación anterior, promediando el desempeño sobre 10 combinaciones de datos distintas y ofreciendo una estimación más robusta y menos dependiente de una división particular de los datos. ## Diagnósticos de residuos

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(modelo_lm, which = 1:4, col = "#3F51B380", pch = 19, cex = 0.8)

par(mfrow = c(1, 1))

Mediante la generación de estos cuatro gráficos de diagnóstico, se evidenció que la mayoría de los supuestos del modelo lineal se cumplen satisfactoriamente. Al analizar el gráfico Residuals vs Fitted, se puede ver que este presenta sus puntos distribuidos aleatoriamente alrededor de cero sin patrones sistemáticos, confirmando el supuesto de linealidad. Luego, el Q-Q plot indica una clara normalidad de residuos, solo presentando desviaciones leves en las colas extremas, lo cual termina por aportar datos insuficientes como para comprometer el supuesto. Y a su vez, la distancia de Cook confirma que ninguna observación ejerce una influencia desproporcionada sobre los coeficientes.

Por último, el único supuesto que podría estar siendo afectado es el de la homocedasticidad, siendo que el gráfico Scale-Location muestra una línea roja con inclinación en ambos extremos, indicando mayor variabilidad en los residuos para tortugas pequeñas y grandes. Aun así, estas desviaciones son poco abruptas y no comprometen la validez general del modelo, por lo que se podría afirmar que los resultados en su totalidad terminaron por cumplir con los supuestos de linealidad planteados.

1.4 Recomendaciones finales

Finalmente, se sugeriría añadir otros predictores que no estén directamente relacionados con la morfología de los individuos, ya que se evidenció que esta posee una gran correlación entre sí y hay un punto en el que deja de brindarnos información relevante para el modelo. Al mismo tiempo, podría ser importante considerar una variable que tenga que ver con el peso de la tortuga o algo de geoestadística, ya que podría alimentar bastante al modelo, considerando que puede que las variaciones morfométricas también estén dadas por el tipo de ambiente en el que nazcan las tortugas. Además de esto, una sugerencia general sería ampliar aún más la base de datos, ya que podría ser bueno para darle al modelo más información y que logre predecir de mejor manera los posibles datos futuros.

1.5 Referencias:

Hermanson, G., & Evers, S. W. (2024). Shell constraints on evolutionary body size-limb size allometry can explain morphological conservatism in the turtle body plan. Ecology and Evolution, 14(11), e70504. https://doi.org/10.1002/ece3.70504 Litzgus, J. D., & Smith, S. E. (2010). Geographic variation in sexual size dimorphism in painted turtles (Chrysemys picta). Journal of Herpetology, 44(2), 320-326. https://doi.org/10.1670/08-333.1