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")

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. En la Parte I no se eliminan; en la Parte II (sección 11) se excluyen para medir su efecto.

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))
)
knitr::kable(comp.tree, digits = 3, caption = "Tabla A. Árbol original vs. árbol podado (MSE en escala log)")
Tabla A. Árbol original vs. árbol podado (MSE en escala log)
Modelo Hojas MSE_train_log MSE_test_log
Árbol original 8 0.136 0.334
Árbol podado (CV) 6 0.157 0.368

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)
knitr::kable(comp.all, digits = 3, caption = "Tabla B. Comparación de modelos: árbol, bagging y bosque aleatorio (MSE de prueba, escala log)")
Tabla B. Comparación de modelos: árbol, bagging y bosque aleatorio (MSE de prueba, escala log)
Modelo MSE_test RMSE_test Mejora_vs_arbol_pct
Árbol (mejor, CV) 0.368 0.607 0.0
Bagging 0.302 0.550 17.8
Bosque Aleatorio 0.285 0.534 22.6

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.

11. Parte II: versión depurada (sin outliers ni variables con ruido)

En la Parte I se usaron todas las observaciones y todas las variables. Aquí se repite el análisis excluyendo lo que genera ruido: (1) jugadores atípicos (como las regiones anómalas y los puntos extremos del scatter de la sección 9) y (2) variables que casi no aportan al bosque. Se usa la misma partición 70/30 y la misma semilla de la Parte I para que la comparación sea justa.

11.1 Detección de outliers (regla 1.5 × IQR)

# (a) Salario atípico dada la trayectoria: residuos de un modelo simple
pre <- lm(logSalary ~ log1p(CAtBat) + log1p(CRuns), data = hitters)
res <- residuals(pre)
out_sal <- res %in% boxplot.stats(res, coef = 1.5)$out

# (b) Trayectorias extremas (superestrellas): CAtBat o CRuns fuera del bigote superior
out_cat <- hitters$CAtBat %in% boxplot.stats(hitters$CAtBat, coef = 1.5)$out
out_cr  <- hitters$CRuns  %in% boxplot.stats(hitters$CRuns,  coef = 1.5)$out

outlier <- out_sal | out_cat | out_cr

tab.out <- data.frame(
  Criterio  = c("Salario atípico dada su trayectoria", "CAtBat extremo", "CRuns extremo", "Total de jugadores excluidos (únicos)"),
  Jugadores = c(sum(out_sal), sum(out_cat), sum(out_cr), sum(outlier))
)
knitr::kable(tab.out, caption = "Tabla 1. Observaciones atípicas detectadas (regla 1.5 × IQR)")
Tabla 1. Observaciones atípicas detectadas (regla 1.5 × IQR)
Criterio Jugadores
Salario atípico dada su trayectoria 5
CAtBat extremo 4
CRuns extremo 8
Total de jugadores excluidos (únicos) 12 
knitr::kable(hitters[outlier, c("CAtBat", "CRuns", "Salary")], caption = "Tabla 2. Jugadores excluidos por ser atípicos (Salary en miles de USD)")
Tabla 2. Jugadores excluidos por ser atípicos (Salary en miles de USD)
CAtBat CRuns Salary
-Don Baylor 7546 1141 950.000
-Darrell Evans 7761 1175 535.000
-Dave Winfield 7287 1135 1861.460
-Graig Nettles 8716 1172 200.000
-Jeffrey Leonard 2964 379 100.000
-Jim Rice 7127 1104 2412.500
-Mike Schmidt 41 6 2127.333
-Pete Rose 14053 2165 750.000
-Reggie Jackson 9528 1509 487.500
-Steve Garvey 8759 1138 1450.000
-Steve Sax 3070 420 90.000
-Terry Kennedy 19 2 920.000 
ggplot(hitters, aes(CAtBat, CRuns, color = outlier)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("gray40", "firebrick"), labels = c("Se conserva", "Outlier excluido")) +
  labs(title = "Outliers en el plano CAtBat vs CRuns", color = NULL) +
  theme_minimal()

11.2 Variables con ruido

Se consideran ruido las variables cuyo %IncMSE en el mejor bosque de la Parte I (entrenado solo con datos de entrenamiento) es menor a 3.

imp.base <- importance(best.forest)[, "%IncMSE"]
ruido <- names(imp.base)[imp.base < 3]

knitr::kable(data.frame(Variable = names(imp.base), IncMSE = round(imp.base, 2),
                        Decision = ifelse(names(imp.base) %in% ruido, "Se excluye", "Se conserva"),
                        row.names = NULL)[order(-imp.base), ],
             caption = "Tabla 3. Variables conservadas y excluidas según su importancia (%IncMSE)")
Tabla 3. Variables conservadas y excluidas según su importancia (%IncMSE)
Variable IncMSE Decision
6 CAtBat 15.82 Se conserva
9 CRuns 14.97 Se conserva
10 CRBI 13.73 Se conserva
7 CHits 13.67 Se conserva
1 AtBat 12.46 Se conserva
11 CWalks 11.94 Se conserva
3 Runs 10.17 Se conserva
8 CHmRun 8.73 Se conserva
5 Walks 8.24 Se conserva
4 RBI 7.12 Se conserva
2 HmRun 4.98 Se conserva
14 PutOuts 2.98 Se excluye
17 NewLeague 2.90 Se excluye
12 League 2.70 Se excluye
16 Errors -0.14 Se excluye
13 Division -0.99 Se excluye
15 Assists -1.33 Se excluye 
dat2    <- dat[, !(names(dat) %in% ruido)]
train2  <- train[!outlier[train]]
test.all   <- setdiff(seq_len(nrow(dat)), train)
test.clean <- test.all[!outlier[test.all]]
keep       <- test.all %in% test.clean

c(entrenamiento = length(train2), prueba_completa = length(test.all), prueba_sin_outliers = length(test.clean))
##       entrenamiento     prueba_completa prueba_sin_outliers 
##                 179                  79                  72

11.3 Árbol, poda y ensambles (versión depurada)

tree2.full <- tree(logSalary ~ ., data = dat2, subset = train2)
set.seed(20261008)
cv2     <- cv.tree(tree2.full)
leaves2 <- sum(tree2.full$frame$var == "<leaf>")
best2   <- cv2$size[which.min(cv2$dev)]
tree2   <- if (best2 < leaves2) prune.tree(tree2.full, best = best2) else tree2.full

plot(cv2$size, cv2$dev, type = "b", xlab = "Número de hojas", ylab = "Deviance (CV)", main = "Validación cruzada (versión depurada)")
points(best2, min(cv2$dev), col = "red", pch = 19, cex = 1.5)

plot(tree2); text(tree2, cex = 0.8, col = "darkgreen")

cat("Hojas del árbol original:", leaves2, "| Tamaño óptimo CV:", best2, "\n")
## Hojas del árbol original: 8 | Tamaño óptimo CV: 8
p2 <- ncol(dat2) - 1
set.seed(20261008)
bag2 <- randomForest(logSalary ~ ., data = dat2, subset = train2, mtry = p2, importance = TRUE)
set.seed(20261008)
rf2  <- randomForest(logSalary ~ ., data = dat2, subset = train2, mtry = max(1, floor(p2 / 3)), importance = TRUE)

yh.tree2 <- predict(tree2, newdata = dat2[test.all, ])
yh.bag2  <- predict(bag2,  newdata = dat2[test.all, ])
yh.rf2   <- predict(rf2,   newdata = dat2[test.all, ])

11.4 Comparación: versión original vs. depurada

y.all <- dat$logSalary[test.all]
mse.pair <- function(yhat) c(completa = mse(yhat, y.all), sin_outliers = mse(yhat[keep], y.all[keep]))

m <- rbind(mse.pair(yhat.best), mse.pair(yhat.bag), mse.pair(yhat.rf),
           mse.pair(yh.tree2),  mse.pair(yh.bag2),  mse.pair(yh.rf2))

comp.v2 <- data.frame(
  Modelo  = rep(c("Árbol podado", "Bagging", "Bosque Aleatorio"), 2),
  Version = rep(c("Original (todas las variables y jugadores)", "Depurada (sin outliers ni ruido)"), each = 3),
  MSE_prueba_completa      = m[, "completa"],
  MSE_prueba_sin_outliers  = m[, "sin_outliers"],
  RMSE_prueba_completa     = sqrt(m[, "completa"])
)
knitr::kable(comp.v2, digits = 3, row.names = FALSE,
             caption = "Tabla 4. MSE de prueba (escala log): modelos originales vs. depurados. 'Completa' = los 79 jugadores de prueba; 'sin outliers' = excluyendo atípicos del set de prueba")
Tabla 4. MSE de prueba (escala log): modelos originales vs. depurados. ‘Completa’ = los 79 jugadores de prueba; ‘sin outliers’ = excluyendo atípicos del set de prueba
Modelo Version MSE_prueba_completa MSE_prueba_sin_outliers RMSE_prueba_completa
Árbol podado Original (todas las variables y jugadores) 0.368 0.224 0.607
Bagging Original (todas las variables y jugadores) 0.302 0.184 0.550
Bosque Aleatorio Original (todas las variables y jugadores) 0.285 0.172 0.534
Árbol podado Depurada (sin outliers ni ruido) 0.359 0.185 0.599
Bagging Depurada (sin outliers ni ruido) 0.365 0.181 0.604
Bosque Aleatorio Depurada (sin outliers ni ruido) 0.348 0.169 0.590

Nota: la columna “completa” es la comparación más exigente, porque ambos modelos se evalúan sobre los mismos jugadores, incluidos los atípicos. La columna “sin outliers” muestra el desempeño en jugadores “típicos”.

11.5 Importancia de variables (versión depurada)

best2.forest <- if (mse(yh.bag2, y.all) <= mse(yh.rf2, y.all)) bag2 else rf2
imp2 <- importance(best2.forest)
imp2.df <- data.frame(Variable = rownames(imp2), imp2, check.names = FALSE)

ggplot(imp2.df, aes(x = reorder(Variable, `%IncMSE`), y = `%IncMSE`)) +
  geom_col(fill = "steelblue") + coord_flip() +
  labs(title = "Importancia de variables (%IncMSE) - versión depurada", x = NULL, y = "% Incremento en MSE") +
  theme_minimal()

knitr::kable(imp2.df[order(-imp2.df$`%IncMSE`), ][1:5, ], digits = 2, row.names = FALSE,
             caption = "Tabla 5. Top 5 de variables más importantes (versión depurada)")
Tabla 5. Top 5 de variables más importantes (versión depurada)
Variable %IncMSE IncNodePurity
CAtBat 15.79 26.19
AtBat 15.33 5.54
CRuns 15.21 22.20
CHits 14.28 20.80
CRBI 12.68 13.19 
top2b <- imp2.df %>% arrange(desc(`%IncMSE`)) %>% slice(1:2) %>% pull(Variable)
top2b
## [1] "CAtBat" "AtBat"

11.6 Árbol de 5 hojas (versión depurada)

f5b <- as.formula(paste("logSalary ~", paste(top2b, collapse = " + ")))
tree.2vb <- tree(f5b, data = dat2, subset = train2,
                 control = tree.control(nobs = length(train2), mindev = 0.001, minsize = 5))
tree.5b <- prune.tree(tree.2vb, best = 5)
summary(tree.5b)
## 
## Regression tree:
## snip.tree(tree = tree.2vb, nodes = c(15L, 4L, 6L, 5L, 14L))
## Number of terminal nodes:  5 
## Residual mean deviance:  0.1321 = 22.98 / 174 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.16200 -0.20780 -0.01385  0.00000  0.22150  1.02000
plot(tree.5b); text(tree.5b, cex = 0.9, col = "darkgreen")

mse.5b <- c(completa = mse(predict(tree.5b, newdata = dat2[test.all, ]), y.all),
            sin_outliers = mse(predict(tree.5b, newdata = dat2[test.clean, ]), dat$logSalary[test.clean]))
knitr::kable(data.frame(Modelo = c("Árbol de 5 hojas (original)", "Árbol de 5 hojas (depurado)"),
                        MSE_prueba_completa = c(mse.tree5, mse.5b["completa"]),
                        MSE_prueba_sin_outliers = c(mse(predict(tree.5, newdata = dat[test.clean, ]), dat$logSalary[test.clean]), mse.5b["sin_outliers"])),
             digits = 3, row.names = FALSE,
             caption = "Tabla 6. Árbol de 5 hojas: versión original vs. depurada (MSE de prueba, escala log)")
Tabla 6. Árbol de 5 hojas: versión original vs. depurada (MSE de prueba, escala log)
Modelo MSE_prueba_completa MSE_prueba_sin_outliers
Árbol de 5 hojas (original) 0.323 0.225
Árbol de 5 hojas (depurado) 0.374 0.229 
vx <- top2b[1]; vy <- top2b[2]
grid <- expand.grid(seq(min(dat2[[vx]]), max(dat2[[vx]]), length.out = 400),
                    seq(min(dat2[[vy]]), max(dat2[[vy]]), length.out = 400))
names(grid) <- c(vx, vy)
grid$pred <- predict(tree.5b, 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)"))
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 = dat2[!outlier, ], 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 del árbol de 5 hojas (versión depurada)", x = vx, y = vy, fill = "Región") +
  theme_minimal()

11.7 Conclusiones de la versión depurada

  • Se excluyeron 12 jugadores atípicos y 6 variables con ruido (League, Division, PutOuts, Assists, Errors, NewLeague).
  • El mejor modelo depurado fue Bosque Aleatorio, con MSE de prueba (completa) de 0.348, contra 0.285 del mejor modelo original.
  • Las dos variables más importantes ahora son CAtBat y AtBat.
  • Excluir outliers puede mejorar el ajuste en jugadores “típicos” pero el modelo pierde capacidad de predecir superestrellas; esta decisión debe justificarse según el objetivo del análisis.