0. Librerías

paquetes <- c("dplyr", "tree", "ISLR", "caret", "ggplot2", "randomForest", "gbm")
faltan   <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltan) > 0) install.packages(faltan, repos = "https://cloud.r-project.org")
## package 'tree' successfully unpacked and MD5 sums checked
## package 'ISLR' successfully unpacked and MD5 sums checked
## package 'gbm' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\52812\AppData\Local\Temp\RtmpmSSdVW\downloaded_packages
library(ISLR)
library(tree)
library(dplyr)
library(ggplot2)
library(caret)
library(randomForest)
library(gbm)

1. Carga y limpieza de datos

El archivo Hitters_in_.csv debe estar en la misma carpeta que este .Rmd. Su primera columna (nombres de los jugadores) viene dañada (#NAME?) y no aporta información, por lo que se elimina. Si el archivo no se encuentra, se usa la base Hitters de la librería ISLR (es la misma base).

if (file.exists("Hitters_in_.csv")) {
  hitters <- read.csv("Hitters_in_.csv", stringsAsFactors = FALSE)
  hitters <- hitters[, -1]            # columna de nombres (#NAME?)
} else {
  data(Hitters, package = "ISLR")
  hitters <- Hitters
}

# Variables categóricas como factor
hitters$League    <- as.factor(hitters$League)
hitters$Division  <- as.factor(hitters$Division)
hitters$NewLeague <- as.factor(hitters$NewLeague)

glimpse(hitters)
## Rows: 322
## Columns: 20
## $ AtBat     <int> 293, 315, 479, 496, 321, 594, 185, 298, 323, 401, 574, 202, …
## $ Hits      <int> 66, 81, 130, 141, 87, 169, 37, 73, 81, 92, 159, 53, 113, 60,…
## $ HmRun     <int> 1, 7, 18, 20, 10, 4, 1, 0, 6, 17, 21, 4, 13, 0, 7, 3, 20, 2,…
## $ Runs      <int> 30, 24, 66, 65, 39, 74, 23, 24, 26, 49, 107, 31, 48, 30, 29,…
## $ RBI       <int> 29, 38, 72, 78, 42, 51, 8, 24, 32, 66, 75, 26, 61, 11, 27, 1…
## $ Walks     <int> 14, 39, 76, 37, 30, 35, 21, 7, 8, 65, 59, 27, 47, 22, 30, 11…
## $ Years     <int> 1, 14, 3, 11, 2, 11, 2, 3, 2, 13, 10, 9, 4, 6, 13, 3, 15, 5,…
## $ CAtBat    <int> 293, 3449, 1624, 5628, 396, 4408, 214, 509, 341, 5206, 4631,…
## $ CHits     <int> 66, 835, 457, 1575, 101, 1133, 42, 108, 86, 1332, 1300, 467,…
## $ CHmRun    <int> 1, 69, 63, 225, 12, 19, 1, 0, 6, 253, 90, 15, 41, 4, 36, 3, …
## $ CRuns     <int> 30, 321, 224, 828, 48, 501, 30, 41, 32, 784, 702, 192, 205, …
## $ CRBI      <int> 29, 414, 266, 838, 46, 336, 9, 37, 34, 890, 504, 186, 204, 1…
## $ CWalks    <int> 14, 375, 263, 354, 33, 194, 24, 12, 8, 866, 488, 161, 203, 2…
## $ League    <fct> A, N, A, N, N, A, N, A, N, A, A, N, N, A, N, A, N, A, A, N, …
## $ Division  <fct> E, W, W, E, E, W, E, W, W, E, E, W, E, E, E, W, W, W, W, W, …
## $ PutOuts   <int> 446, 632, 880, 200, 805, 282, 76, 121, 143, 0, 238, 304, 211…
## $ Assists   <int> 33, 43, 82, 11, 40, 421, 127, 283, 290, 0, 445, 45, 11, 151,…
## $ Errors    <int> 20, 10, 14, 3, 4, 25, 7, 9, 19, 0, 22, 11, 7, 6, 8, 0, 10, 1…
## $ Salary    <dbl> NA, 475.000, 480.000, 500.000, 91.500, 750.000, 70.000, 100.…
## $ NewLeague <fct> A, N, A, N, N, A, A, A, N, A, A, N, N, A, N, A, N, A, A, N, …

Según las instrucciones de la actividad, se excluyen las variables Years y Hits.

hitters <- hitters %>% dplyr::select(-Years, -Hits)

cat("Filas duplicadas:", sum(duplicated(hitters)), "\n")
## Filas duplicadas: 0
colSums(is.na(hitters))
##     AtBat     HmRun      Runs       RBI     Walks    CAtBat     CHits    CHmRun 
##         0         0         0         0         0         0         0         0 
##     CRuns      CRBI    CWalks    League  Division   PutOuts   Assists    Errors 
##         0         0         0         0         0         0         0         0 
##    Salary NewLeague 
##        59         0
summary(hitters$Salary)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    67.5   190.0   425.0   535.9   750.0  2460.0      59

Los únicos valores faltantes están en la variable dependiente (Salary). Como no se puede entrenar ni evaluar un modelo sin la variable objetivo, se eliminan esas observaciones.

cat("Filas antes:", nrow(hitters), "\n")
## Filas antes: 322
hitters <- na.omit(hitters)
cat("Filas después:", nrow(hitters), "\n")
## Filas después: 263
colSums(is.na(hitters))
##     AtBat     HmRun      Runs       RBI     Walks    CAtBat     CHits    CHmRun 
##         0         0         0         0         0         0         0         0 
##     CRuns      CRBI    CWalks    League  Division   PutOuts   Assists    Errors 
##         0         0         0         0         0         0         0         0 
##    Salary NewLeague 
##         0         0

2. Análisis Exploratorio de Datos

2.1 Histogramas

par(mfrow = c(1, 2))
hist(hitters$Salary, main = "Salary (miles de USD)", xlab = "Salary", col = "lightsteelblue")
hist(log(hitters$Salary), main = "log(Salary)", xlab = "log(Salary)", col = "darkseagreen")

par(mfrow = c(1, 1))

Salary tiene un sesgo positivo marcado (muchos salarios bajos y pocos muy altos). Con el logaritmo la distribución se parece mucho más a una normal, por lo que se transforma la variable objetivo a log(Salary).

hitters$logSalary <- log(hitters$Salary)
num_pred <- hitters %>% dplyr::select(where(is.numeric), -Salary, -logSalary)

par(mfrow = c(4, 4), mar = c(3, 3, 2, 1))
for (v in names(num_pred)) hist(num_pred[[v]], main = v, xlab = "", col = "lightsteelblue")
par(mfrow = c(1, 1))

Las estadísticas de carrera (CAtBat, CHits, CRuns, etc.) también son asimétricas a la derecha. Esto no requiere transformarlas, porque los árboles de decisión solo usan el orden de los valores y son invariantes a transformaciones monótonas de las variables independientes.

2.2 Matriz de correlaciones

num_all <- hitters %>% dplyr::select(where(is.numeric))
cm <- cor(num_all)

cm_df <- as.data.frame(as.table(cm))
ggplot(cm_df, aes(Var1, Var2, fill = Freq)) +
  geom_tile() +
  geom_text(aes(label = round(Freq, 2)), size = 2.3) +
  scale_fill_gradient2(low = "firebrick", mid = "white", high = "steelblue", limits = c(-1, 1)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(x = NULL, y = NULL, fill = "r", title = "Matriz de correlaciones")

sort(cm[, "logSalary"], decreasing = TRUE)
##   logSalary      Salary       CRuns       CHits      CAtBat        CRBI 
##  1.00000000  0.89543720  0.62108935  0.61994113  0.61079998  0.60267128 
##      CWalks      CHmRun         RBI       Walks        Runs       AtBat 
##  0.54658066  0.52394269  0.44415019  0.43243922  0.42559551  0.41494587 
##       HmRun     PutOuts     Assists      Errors 
##  0.33985428  0.22448262  0.04996914 -0.02076464

2.3 Boxplots

par(mfrow = c(1, 3))
boxplot(logSalary ~ League,    data = hitters, col = "lightsteelblue", main = "log(Salary) por League")
boxplot(logSalary ~ Division,  data = hitters, col = "lightsteelblue", main = "log(Salary) por Division")
boxplot(logSalary ~ NewLeague, data = hitters, col = "lightsteelblue", main = "log(Salary) por NewLeague")

par(mfrow = c(1, 1))
# Variables estandarizadas para poder compararlas en una misma escala
boxplot(scale(num_pred), las = 2, col = "darkseagreen", main = "Variables numéricas (estandarizadas)")

Se observan varios valores atípicos en las estadísticas de carrera. No se eliminan, ya que corresponden a jugadores veteranos reales y los árboles son robustos a ellos.

3. Separación en entrenamiento y prueba (70/30)

dat <- hitters %>% dplyr::select(-Salary)      # logSalary será la variable objetivo

set.seed(20261008)
train <- sample(1:nrow(dat), round(nrow(dat) * 0.7))

dat.test    <- dat[-train, ]
y.test      <- dat$logSalary[-train]
salary.test <- hitters$Salary[-train]

c(entrenamiento = length(train), prueba = nrow(dat.test))
## entrenamiento        prueba 
##           184            79

4. Árbol de decisión con todas las variables

tree.full <- tree(logSalary ~ ., data = dat, subset = train)
summary(tree.full)
## 
## Regression tree:
## tree(formula = logSalary ~ ., data = dat, subset = train)
## Variables actually used in tree construction:
## [1] "CAtBat" "CRuns"  "AtBat"  "CWalks" "CRBI"   "Walks" 
## Number of terminal nodes:  8 
## Residual mean deviance:  0.1423 = 25.05 / 176 
## Distribution of residuals:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.558000 -0.193700 -0.006555  0.000000  0.230200  1.485000
plot(tree.full)
text(tree.full, cex = 0.75, col = "darkgreen")

tree.full
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 184 138.4000 5.937  
##    2) CAtBat < 1452 72  21.1400 5.099  
##      4) CRuns < 91 37   7.6630 4.756  
##        8) AtBat < 158 5   3.0500 5.339 *
##        9) AtBat > 158 32   2.6460 4.665 *
##      5) CRuns > 91 35   4.5270 5.461  
##       10) CWalks < 116 24   1.7440 5.301 *
##       11) CWalks > 116 11   0.8247 5.810 *
##    3) CAtBat > 1452 112  34.0100 6.477  
##      6) AtBat < 421 46   7.1110 6.149 *
##      7) AtBat > 421 66  18.5300 6.705  
##       14) CRBI < 273 12   3.3770 6.058 *
##       15) CRBI > 273 54   9.0030 6.849  
##         30) Walks < 65.5 38   4.4170 6.704 *
##         31) Walks > 65.5 16   1.8830 7.194 *

5. Validación cruzada y poda

set.seed(20261008)
cv.hit <- cv.tree(tree.full)
cv.hit
## $size
## [1] 8 7 6 5 4 3 2 1
## 
## $dev
## [1]  46.11128  44.24910  43.94693  45.27242  53.89075  64.03955  66.30191
## [8] 141.21159
## 
## $k
## [1]      -Inf  1.958075  1.967368  2.702826  6.145657  8.372023  8.945613
## [8] 83.230454
## 
## $method
## [1] "deviance"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"
plot(cv.hit$size, cv.hit$dev, type = "b", xlab = "Número de hojas", ylab = "Deviance (CV)",
     main = "Validación cruzada")

n.leaves  <- sum(tree.full$frame$var == "<leaf>")
best.size <- cv.hit$size[which.min(cv.hit$dev)]
points(best.size, min(cv.hit$dev), col = "red", pch = 19, cex = 1.5)

cat("Hojas del árbol original:", n.leaves, "\nTamaño óptimo según CV:", best.size, "\n")
## Hojas del árbol original: 8 
## Tamaño óptimo según CV: 6
# Si la CV sugiere menos hojas, se poda; si no, el árbol original ya es el mejor
tree.best <- if (best.size < n.leaves) prune.tree(tree.full, best = best.size) else tree.full
plot(tree.best)
text(tree.best, cex = 0.8, col = "darkgreen")

Comparación: árbol original vs. árbol podado

mse <- function(pred, real) mean((pred - real)^2)

yhat.full <- predict(tree.full, newdata = dat.test)
yhat.best <- predict(tree.best, newdata = dat.test)

comp.tree <- data.frame(
  Modelo        = c("Árbol original", "Árbol podado (CV)"),
  Hojas         = c(n.leaves, sum(tree.best$frame$var == "<leaf>")),
  MSE_train_log = c(mse(predict(tree.full, dat[train, ]), dat$logSalary[train]),
                    mse(predict(tree.best, dat[train, ]), dat$logSalary[train])),
  MSE_test_log  = c(mse(yhat.full, y.test), mse(yhat.best, y.test))
)
comp.tree

6. Predicciones en datos de prueba (mejor árbol)

mse.tree <- mse(yhat.best, y.test)

ggplot(data.frame(real = y.test, pred = yhat.best), aes(x = pred, y = real)) +
  geom_point(color = "steelblue", size = 2) +
  geom_abline(intercept = 0, slope = 1, color = "red") +
  labs(title = "Árbol de decisión: predicho vs. real (prueba)",
       x = "log(Salary) predicho", y = "log(Salary) real") +
  theme_minimal()

cat("MSE de prueba (escala log):", round(mse.tree, 4), "\n")
## MSE de prueba (escala log): 0.3679
cat("RMSE (escala log):", round(sqrt(mse.tree), 4), "\n")
## RMSE (escala log): 0.6065
# En dólares (miles de USD), regresando con exp()
rmse.tree.usd <- sqrt(mse(exp(yhat.best), salary.test))
cat("RMSE en escala original (miles de USD):", round(rmse.tree.usd, 1), "\n")
## RMSE en escala original (miles de USD): 392.4

El error principal se reporta en escala logarítmica porque es la variable con la que se entrenó el modelo. El RMSE en miles de USD es solo una referencia aproximada, ya que exp() de una predicción en logaritmos estima la mediana y no la media del salario.

7. Bagging y Bosques Aleatorios

p <- ncol(dat) - 1        # número de variables independientes

# Bagging: mtry = p (todas las variables en cada partición)
set.seed(20261008)
bag.hit <- randomForest(logSalary ~ ., data = dat, subset = train, mtry = p, importance = TRUE)
bag.hit
## 
## Call:
##  randomForest(formula = logSalary ~ ., data = dat, mtry = p, importance = TRUE,      subset = train) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 17
## 
##           Mean of squared residuals: 0.1798356
##                     % Var explained: 76.09
# Bosque Aleatorio: mtry = p/3 (valor por default en regresión)
set.seed(20261008)
rf.hit <- randomForest(logSalary ~ ., data = dat, subset = train, mtry = floor(p / 3), importance = TRUE)
rf.hit
## 
## Call:
##  randomForest(formula = logSalary ~ ., data = dat, mtry = floor(p/3),      importance = TRUE, subset = train) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##           Mean of squared residuals: 0.1699338
##                     % Var explained: 77.4
yhat.bag <- predict(bag.hit, newdata = dat.test)
yhat.rf  <- predict(rf.hit,  newdata = dat.test)

mse.bag <- mse(yhat.bag, y.test)
mse.rf  <- mse(yhat.rf,  y.test)

par(mfrow = c(1, 2))
plot(yhat.bag, y.test, main = "Bagging", xlab = "Predicho", ylab = "Real", pch = 19, col = "steelblue"); abline(0, 1, col = "red")
plot(yhat.rf,  y.test, main = "Bosque Aleatorio", xlab = "Predicho", ylab = "Real", pch = 19, col = "darkseagreen"); abline(0, 1, col = "red")

par(mfrow = c(1, 1))

Comparación de modelos

comp.all <- data.frame(
  Modelo   = c("Árbol (mejor, CV)", "Bagging", "Bosque Aleatorio"),
  MSE_test = c(mse.tree, mse.bag, mse.rf)
)
comp.all$RMSE_test <- sqrt(comp.all$MSE_test)
comp.all$Mejora_vs_arbol_pct <- round(100 * (1 - comp.all$MSE_test / mse.tree), 1)
comp.all

8. Importancia de las variables

best.forest <- if (mse.bag <= mse.rf) bag.hit else rf.hit
best.name   <- if (mse.bag <= mse.rf) "Bagging" else "Bosque Aleatorio"
cat("Modelo con menor MSE de prueba:", best.name, "\n")
## Modelo con menor MSE de prueba: Bosque Aleatorio
imp <- importance(best.forest)
imp
##              %IncMSE IncNodePurity
## AtBat     12.4590573     6.2368488
## HmRun      4.9831730     1.8980108
## Runs      10.1660902     3.7537964
## RBI        7.1202754     3.9866506
## Walks      8.2436397     4.9799545
## CAtBat    15.8160976    25.7529094
## CHits     13.6669573    20.5504346
## CHmRun     8.7320880     6.9747041
## CRuns     14.9735922    24.4318064
## CRBI      13.7303875    15.2475814
## CWalks    11.9393879    17.1717467
## League     2.6993311     0.3432669
## Division  -0.9853962     0.2458668
## PutOuts    2.9819112     1.7560706
## Assists   -1.3278166     1.2157146
## Errors    -0.1352357     1.1280919
## NewLeague  2.8970751     0.4427656
imp.df <- data.frame(Variable = rownames(imp), imp, check.names = FALSE)

ggplot(imp.df, aes(x = reorder(Variable, `%IncMSE`), y = `%IncMSE`)) +
  geom_col(fill = "steelblue") + coord_flip() +
  labs(title = paste("Importancia de variables (%IncMSE) -", best.name), x = NULL, y = "% Incremento en MSE") +
  theme_minimal()

ggplot(imp.df, aes(x = reorder(Variable, IncNodePurity), y = IncNodePurity)) +
  geom_col(fill = "darkseagreen") + coord_flip() +
  labs(title = paste("Importancia de variables (IncNodePurity) -", best.name), x = NULL, y = "Incremento en pureza de nodos") +
  theme_minimal()

varImpPlot(best.forest, main = paste("varImpPlot -", best.name))

top2 <- imp.df %>% arrange(desc(`%IncMSE`)) %>% slice(1:2) %>% pull(Variable)
top2
## [1] "CAtBat" "CRuns"

9. Árbol de 5 hojas con las 2 variables más relevantes

f5 <- as.formula(paste("logSalary ~", paste(top2, collapse = " + ")))

# Se deja crecer el árbol lo suficiente para poder podarlo a exactamente 5 hojas
tree.2v <- tree(f5, data = dat, subset = train,
                control = tree.control(nobs = length(train), mindev = 0.001, minsize = 5))
tree.5  <- prune.tree(tree.2v, best = 5)

summary(tree.5)
## 
## Regression tree:
## snip.tree(tree = tree.2v, nodes = c(5L, 9L, 7L, 6L))
## Number of terminal nodes:  5 
## Residual mean deviance:  0.2088 = 37.38 / 179 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.86000 -0.29300  0.01791  0.00000  0.29360  1.21100
tree.5
## node), split, n, deviance, yval
##       * denotes terminal node
## 
## 1) root 184 138.400 5.937  
##   2) CAtBat < 1452 72  21.140 5.099  
##     4) CRuns < 91 37   7.663 4.756  
##       8) CRuns < 17 2   1.244 6.036 *
##       9) CRuns > 17 35   2.955 4.683 *
##     5) CRuns > 91 35   4.527 5.461 *
##   3) CAtBat > 1452 112  34.010 6.477  
##     6) CRuns < 763 87  22.160 6.359 *
##     7) CRuns > 763 25   6.492 6.884 *
plot(tree.5)
text(tree.5, cex = 0.9, col = "darkgreen")

mse.tree5 <- mse(predict(tree.5, newdata = dat.test), y.test)
cat("MSE de prueba del árbol de 5 hojas (escala log):", round(mse.tree5, 4), "\n")
## MSE de prueba del árbol de 5 hojas (escala log): 0.3226

Scatter plot con las 5 regiones

Las regiones se obtienen evaluando el árbol de 5 hojas sobre una malla fina del plano formado por las dos variables; cada región es una hoja del árbol y se etiqueta con el salario medio que predice (en miles de USD).

vx <- top2[1]; vy <- top2[2]

grid <- expand.grid(seq(min(dat[[vx]]), max(dat[[vx]]), length.out = 400),
                    seq(min(dat[[vy]]), max(dat[[vy]]), length.out = 400))
names(grid) <- c(vx, vy)
grid$pred <- predict(tree.5, newdata = grid)

lv <- sort(unique(round(grid$pred, 4)))
grid$region <- factor(round(grid$pred, 4), levels = lv,
                      labels = paste0("Región ", seq_along(lv), " (~$", round(exp(lv)), "k)"))

# Posición de la etiqueta de cada región
lab <- grid %>% group_by(region) %>%
  summarise(x = mean(.data[[vx]]), y = mean(.data[[vy]]), .groups = "drop")
lab$id <- paste0("R", seq_len(nrow(lab)))

ggplot() +
  geom_raster(data = grid, aes(x = .data[[vx]], y = .data[[vy]], fill = region), alpha = 0.40) +
  geom_point(data = dat, aes(x = .data[[vx]], y = .data[[vy]]), size = 1.8, color = "gray20") +
  geom_text(data = lab, aes(x = x, y = y, label = id), fontface = "bold", size = 6) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Regiones (hojas) del árbol de 5 hojas", x = vx, y = vy, fill = "Región") +
  theme_minimal()

10. Conclusiones

  • Limpieza y transformación. Se eliminaron las observaciones sin Salary y se trabajó con log(Salary) para corregir el sesgo de la variable objetivo.
  • Árbol de decisión. La validación cruzada sugirió un árbol de 6 hojas, contra 8 del árbol original. Su MSE de prueba (escala log) fue de 0.368.
  • Modelos de ensamble. El Bagging obtuvo un MSE de prueba de 0.302 y el Bosque Aleatorio de 0.285. El mejor modelo en prueba fue: Bosque Aleatorio. Combinar muchos árboles redujo el error respecto a un solo árbol, ya que promediar árboles disminuye la varianza que tiene un árbol individual.
  • Variables más relevantes. Las dos variables más importantes fueron CAtBat y CRuns.
  • Árbol de 5 hojas. Con solo esas dos variables se obtuvo un modelo mucho más sencillo de interpretar, con MSE de prueba de 0.323, comparable contra los modelos anteriores y con regiones fáciles de visualizar en el scatter plot.
  • Trade-off. Un solo árbol es muy interpretable pero más inestable; Bagging y Bosques Aleatorios predicen mejor pero se pierde la lectura directa de las reglas. Para explicar y comunicar sirve el árbol de 5 hojas; para predecir, el bosque.

Nota: revisa las gráficas y tablas ya con los resultados y ajusta el texto de interpretación (por ejemplo, qué significan en términos de béisbol las variables más importantes y qué tipo de jugador cae en cada región) antes de presentarlo.