Regresión logística

En la clase anterior trabajamos con regresión lineal, donde la variable respuesta era un número continuo (ventas). Cerramos esa clase con una pregunta:

¿Qué pasa cuando la variable que queremos predecir no es un número, sino una categoría?

Hoy vamos a responder esa pregunta con el modelo más simple para este tipo de problemas: la regresión logística.

El problema

El dataset Ingresos contiene información de 500 clientes: su ingreso mensual (en miles) y si terminaron adquiriendo o no un producto.

Nuestra pregunta:

¿Podemos predecir si un cliente va a adquirir el producto a partir de su ingreso mensual?

library(readxl)
library(ggplot2)

ingresos <- read_excel("Ingresos.xlsx")

head(ingresos)
## # A tibble: 6 × 2
##   Ingresos_mensuales_k Adquiere_producto
##                  <dbl>             <dbl>
## 1                 21.9                 0
## 2                 47.8                 1
## 3                 37.9                 1
## 4                 31.9                 1
## 5                 12.0                 0
## 6                 12.0                 0


En este caso la variable respuesta (Adquiere_producto) ya no es continua: solo toma dos valores, 0 (no adquiere) o 1 (sí adquiere). Le decimos a R explícitamente que la trate como una categoría y no como un número:

ingresos$Adquiere_producto <- as.factor(ingresos$Adquiere_producto)

str(ingresos)
## tibble [500 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Ingresos_mensuales_k: num [1:500] 21.9 47.8 37.9 31.9 12 ...
##  $ Adquiere_producto   : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 2 2 1 ...


Exploración

Como siempre, antes de modelar miramos los datos. Comparemos la distribución del ingreso mensual entre quienes adquirieron el producto y quienes no.

# ggplot() inicia la construcción del gráfico.
#
# Primer argumento: el data.frame que contiene los datos.
# aes() (aesthetics) define cómo se asignan las variables a los elementos
# visuales del gráfico.
#
# En este caso:
# - x: la variable que irá sobre el eje horizontal.
# - y: la variable que irá sobre el eje vertical.
# - fill: el color de relleno de las cajas, que también dependerá de
#         la variable Adquiere_producto.
ggplot(
  ingresos,
  aes(
    x = Adquiere_producto,
    y = Ingresos_mensuales_k,
    fill = Adquiere_producto
  )
) +

  # Primera capa: dibuja un boxplot para cada categoría del eje X.
  geom_boxplot() +

  # Segunda capa: agrega etiquetas descriptivas al gráfico.
  #
  # title: título principal.
  # x: etiqueta del eje X.
  # y: etiqueta del eje Y.
  labs(
    title = "Ingreso mensual según si adquiere el producto",
    x = "¿Adquiere el producto?",
    y = "Ingreso mensual (miles)"
  ) +

  # Tercera capa: modifica la apariencia general del gráfico.
  #
  # Como el color (fill) representa exactamente la misma variable que
  # aparece en el eje X, la leyenda sería redundante y se elimina.
  theme(
    legend.position = "none"
  )


Se nota que quienes adquieren el producto tienden a tener ingresos más altos. Eso ya nos dice que la variable Ingresos_mensuales_k probablemente sirva para predecir Adquiere_producto. Con esa señal visual nos alcanza para avanzar al modelo, en este curso no vamos a emplear un test estadístico formal para confirmarlo, como podría ser la prueba de razón de verosimilitud (likelihood test ratio), que nos sirve para evaluar si incluir Ingresos_mensuales_k mejora significativamente el modelo respecto al modelo nulo.

# Modelo nulo
modelo_nulo <- glm(Adquiere_producto ~ 1,
                   data = ingresos,
                   family = binomial)

# Modelo con la variable explicativa
modelo_logit <- glm(Adquiere_producto ~ Ingresos_mensuales_k,
                    data = ingresos,
                    family = binomial)

# Likelihood Ratio Test
anova(modelo_nulo, modelo_logit, test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: Adquiere_producto ~ 1
## Model 2: Adquiere_producto ~ Ingresos_mensuales_k
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       499     666.93                          
## 2       498     414.94  1   251.99 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#  Pr(>Chi) bajo nos dice que sí mejora


¿Por qué no usar lm()?

Podríamos preguntarnos, ¿por qué no usamos regresión lineal, poniendo 0 y 1 como si fueran números? Probemos y veamos qué pasa.

# Crear una copia del conjunto de datos original
ingresos_num <- ingresos

# Convertir la variable "Adquiere_producto" a formato numérico.
# Primero se transforma a carácter para evitar problemas si originalmente es un factor.
ingresos_num$Adquiere_producto <- as.numeric(as.character(ingresos$Adquiere_producto))

# Ajustar un modelo de regresión lineal simple.
# La variable dependiente es "Adquiere_producto" y la variable explicativa es
# "Ingresos_mensuales_k", utilizando el conjunto de datos "ingresos_num".
modelo_lineal <- lm(Adquiere_producto ~ Ingresos_mensuales_k, data = ingresos_num)

# Calcular los valores predichos por el modelo y mostrar su rango
# (valor mínimo y máximo de las predicciones).
range(predict(modelo_lineal))
## [1] -0.1359513  0.9089063


El modelo lineal predice valores como -0.14 o 0.91 que no representan una probabilidad válida (una probabilidad siempre está entre 0 y 1). Para un ingreso muy bajo, incluso podría llegar a predecir un valor negativo. Esto no tiene sentido si lo que queremos es interpretar el resultado como “la probabilidad de que el cliente adquiera el producto”.

Necesitamos una función que, sin importar qué valores tome X, siempre devuelva un número entre 0 y 1. Ahí es donde entra la regresión logística.

La función logística

La regresión logística no ajusta una recta; ajusta una curva en forma de “S” (llamada función logística o sigmoide) que está diseñada específicamente para mantenerse siempre entre 0 y 1:

\[ p(X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X)}} \]

No necesitamos memorizar esta fórmula ni manipularla a mano porque R se encarga de estimar los coeficientes β0 y β1. Lo importante es la forma de la curva:

curve(1 / (1 + exp(-(-5 + 0.15 * x))),
      from = 0, to = 60,
      lwd = 2.5, col = "firebrick",
      xlab = "X (variable predictora)",
      ylab = "Probabilidad estimada",
      main = "Forma de la función logística")
abline(h = c(0,1), lty = 2, col = "black")


Sin importar qué tan extremo sea el valor de X, la curva nunca baja de 0 ni sube de 1. Eso es exactamente lo que necesitábamos.

Ajustando el modelo con glm()

En R, la regresión logística se ajusta con glm() (generalized linear model), agregando el argumento family = binomial para indicarle que la variable respuesta es de dos categorías.

Nota: si la variable respuesta presenta más de dos categorías, es necesario utilizar otras variantes de la regresión logística, como logística multinomial cuando las categorías no tienen un orden natural (e.g., tipo de transporte: auto, autobús o bicicleta) o regresión logística ordinal cuando las categorías poseen un orden (e.g., nivel de satisfacción: bajo, medio y alto).

modelo <- glm(Adquiere_producto ~ Ingresos_mensuales_k, data = ingresos, family = binomial)

summary(modelo)
## 
## Call:
## glm(formula = Adquiere_producto ~ Ingresos_mensuales_k, family = binomial, 
##     data = ingresos)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -4.96764    0.43054  -11.54   <2e-16 ***
## Ingresos_mensuales_k  0.14949    0.01285   11.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 666.93  on 499  degrees of freedom
## Residual deviance: 414.94  on 498  degrees of freedom
## AIC: 418.94
## 
## Number of Fisher Scoring iterations: 5


Igual que con lm(), nos interesan sobre todo dos cosas de esta salida: el signo del coeficiente y la columna Pr(>|z|).

  • El coeficiente de Ingresos_mensuales_k es positivo. Esto nos dice que a mayor ingreso, mayor es la probabilidad estimada de adquirir el producto.
  • Al igual que en regresión lineal, mientras más chico el valor de Pr(>|z|), más confianza tenemos en que esa variable realmente aporta información.

Nota: el coeficiente de la regresión logística no se interpreta directamente como “cuánto sube la probabilidad”, como sí hacíamos con la pendiente en regresión lineal. Se interpreta en una escala distinta (el logit).

Si los coeficientes de la regresión logística están expresados en términos de log-odds, también llamados logit, entonces repasemos… ¿qué escalas podemos usar?

  • Probabilidad

Lo interpretamos como la probabilidad de que ocurra el evento:

\[ p=P(Y=1) \]

Por ejemplo, si la probabilidad de adquirir el producto es 50 %:

\[ p=0.50 \]


  • Odds

Los odds representan la razón entre la probabilidad de que ocurra un evento y la probabilidad de que no ocurra:

\[ \text{Odds} = \frac{P(Y=1)}{1 - P(Y=1)} \]

Si \(p=0.50\):

\[ \text{Odds}=\frac{0.50}{1-0.50}=1 \] Un valor mayor que 1 indica que los momios aumentan, un valor menor que 1 indica que disminuyen y un valor igual a 1 significa que no cambian.


  • Log-odds o logit

La regresión logística trabaja con el logaritmo natural de los odds,

\[ \ln\left(\frac{P(Y=1)}{1-P(Y=1)}\right) \]

Por ejemplo, si \(p=0.50\):

\[ \ln\left(\frac{0.50}{1-0.50}\right) = \ln(1) = 0 \]

Así,

\[ P(Y=1)=0.50 \quad\longrightarrow\quad \text{Odds}=1 \quad\longrightarrow\quad \text{Log-odds}=0 \] Estas tres cantidades representan la misma situación, pero en diferentes escalas.

Y la transformación también funciona en sentido contrario.

Si conocemos el logit, podemos recuperar primero los odds:

\[ \text{Odds}=e^{\text{logit}} \] y después convertir los odds en probabilidad:

\[ P(Y=1)=\frac{\text{Odds}}{1+\text{Odds}} \] Resultando:

\[ P(Y=1)= \frac{e^{\text{logit}}} {1+e^{\text{logit}}} \] Siendo esta es la función logística que produce la conocida curva en forma de S.

Una forma práctica para interpretarlos es a través del cálculo de exp(coeficiente).

exp(coef(modelo)["Ingresos_mensuales_k"])
## Ingresos_mensuales_k 
##             1.161245


Esto significa que por cada mil unidades monetarias adicionales de ingreso mensual, los momios de adquirir el producto se multiplican por 1.16, es decir, aumentan aproximadamente un 16 %. Esto no significa que la probabilidad de adquirir el producto aumente un 16 %, ya que la relación entre probabilidad y momios no es lineal.

Por ejemplo, si inicialmente la probabilidad de adquirir el producto fuera del 50 %, los momios serían:

\[ \text{Odds} = \frac{0.5}{1-0.5} = 1 \] Al multiplicarlos por 1.16, los nuevos momios serían 1.16.

\[ \text{Odds}_{\text{nuevo}} = 1 \times 1.16 = 1.16 \] La probabilidad correspondiente pasa a ser: \[ p_{\text{nuevo}} = \frac{\text{Odds}_{\text{nuevo}}} {1+\text{Odds}_{\text{nuevo}}} = \frac{1.16}{1+1.16} \approx 0.537 \]

Es decir, la probabilidad aumenta de 50 % a 53.7 %, un incremento de 3.7 puntos porcentuales, y no de un 16 %. El cambio exacto en la probabilidad siempre depende del valor inicial de esta.

Visualizando el modelo ajustado

# ggplot() crea la base del gráfico.
#
# Se especifica el data.frame (ingresos) y las variables que se asignarán
# a los ejes mediante aes() (aesthetics).
#
# - x: ingreso mensual.
# - y: la variable respuesta. Como es un factor ("0" y "1"), se convierte
#      a valores numéricos para poder representarla en el eje Y.
ggplot(
  ingresos,
  aes(
    x = Ingresos_mensuales_k,
    y = as.numeric(as.character(Adquiere_producto))
  )
) +

  # Primera capa: dibuja un punto por cada observación.
  #
  # alpha = 0.5 hace los puntos semitransparentes para que sea más fácil
  # apreciar zonas donde muchos puntos se superponen.
  geom_point(
    color = "darkblue",
    alpha = 0.5
  ) +

  # Segunda capa: ajusta un modelo de regresión logística y dibuja
  # la curva de probabilidades estimadas.
  #
  # method = "glm" indica que se utilizará un Modelo Lineal Generalizado.
  # family = "binomial" especifica que se trata de una regresión logística.
  # se = FALSE evita mostrar la banda de confianza alrededor de la curva.
  geom_smooth(
    method = "glm",
    method.args = list(family = "binomial"),
    se = FALSE,
    color = "firebrick"
  ) +

  # Tercera capa: agrega el título y las etiquetas de los ejes.
  labs(
    title = "Probabilidad de adquirir el producto según ingreso mensual",
    x = "Ingreso mensual (miles)",
    y = "P(Adquiere producto = 1)"
  )
## `geom_smooth()` using formula = 'y ~ x'


Ahí está la curva en forma de S que mencionamos antes, ahora ajustada a nuestros datos reales.

Prediciendo probabilidades

Nuestro modelo es: \[ -4.96764 + 0.14949(\text{Ingresos_mensuales_k}) \] Supongamos que queremos estimar la probabilidad de adquirir el producto para una persona con un ingreso mensual de 25 mil unidades monetarias.

-4.96764+0.14949*25 # este valor es el logit
## [1] -1.23039
predict(modelo, newdata = data.frame(Ingresos_mensuales_k = 25))
##         1 
## -1.230317


Para entender qué significa ese -1.23, podemos regresar de la escala de log-odds a la escala de odds:

\[ \text{Odds}=e^{-1.23} \approx 0.292 \] Y convertir los odds en probabilidad:

\[ P(Y=1)= \frac{0.292}{1+0.292} \approx 0.226 \] Por lo tanto, para un ingreso mensual de 25 mil, el modelo estima una probabilidad de aproximadamente 22.6 %.

En R podemos hacer toda esta transformación directamente, con type = "response" le pedimos a R que nos devuelva directamente la probabilidad (no el logit).

# type = "link"      → logit / log-odds
# type = "response"  → probabilidad

predict(modelo, newdata = data.frame(Ingresos_mensuales_k = 25), type = "response")
##        1 
## 0.226126
predict(modelo, newdata = data.frame(Ingresos_mensuales_k = 45), type = "response")
##         1 
## 0.8531526


El logit puede tomar cualquier valor entre \(-\infty\) e \(\infty\), mientras que la probabilidad siempre está entre 0 y 1.

De probabilidad a clasificación

El modelo nos da una probabilidad, pero muchas veces necesitamos una decisión binaria: ¿le ofrezco el producto a este cliente o no? Para eso definimos un umbral (el más simple y común es 0.5): si la probabilidad predicha supera el umbral, clasificamos como 1; si no, como 0.

probabilidad <- predict(modelo, newdata = data.frame(Ingresos_mensuales_k = 25), type = "response")

unname(ifelse(probabilidad >= 0.5, 1, 0))
## [1] 0
# unname simplemente elimina cualquier nombre que pudiera tener ese resultado


Train / test, otra vez

Igual que en regresión lineal, para evaluar el modelo de forma honesta necesitamos separar datos de entrenamiento y de prueba.

set.seed(123)

n <- nrow(ingresos)
indices_entrenamiento <- sample(1:n, size = 0.8 * n)

train <- ingresos[indices_entrenamiento, ]
test <- ingresos[-indices_entrenamiento, ]

modelo_train <- glm(Adquiere_producto ~ Ingresos_mensuales_k, data = train, family = binomial)


Predecimos sobre test (datos que el modelo nunca vio) y clasificamos con el umbral de 0.5:

prob_test <- predict(modelo_train, newdata = test, type = "response")
pred_clase <- ifelse(prob_test >= 0.5, 1, 0)

head(data.frame(real = test$Adquiere_producto, probabilidad = round(prob_test, 3), prediccion = pred_clase),10)
##    real probabilidad prediccion
## 1     0        0.165          0
## 2     0        0.052          0
## 3     0        0.110          0
## 4     1        0.224          0
## 5     0        0.334          0
## 6     0        0.110          0
## 7     0        0.030          0
## 8     1        0.875          1
## 9     1        0.574          1
## 10    0        0.115          0


Matriz de confusión

Para regresión lineal usábamos RMSE para medir el error. Para clasificación, el punto de partida es la matriz de confusión: una tabla que cruza lo que el modelo predijo contra lo que realmente pasó.

La matriz de confusión es una herramienta que permite evaluar el desempeño de un modelo de clasificación comparando las predicciones del modelo con los valores reales. En ella se organizan los resultados en una tabla donde las filas representan los valores reales y las columnas las predicciones del modelo.

matriz_confusion <- table(Real = test$Adquiere_producto, Predicho = pred_clase)
matriz_confusion
##     Predicho
## Real  0  1
##    0 52 12
##    1  6 30


library(ggplot2)

# Convertir la matriz de confusión en un data frame para poder utilizar ggplot2
df <- as.data.frame(matriz_confusion)

# Definir el orden de las categorías de las variables
# (0 antes que 1 tanto para los valores reales como para las predicciones)
df$Real <- factor(df$Real, levels = c("0", "1"))
df$Predicho <- factor(df$Predicho, levels = c("0", "1"))

# Representar la matriz de confusión como un mapa de calor
ggplot(df, aes(Predicho, Real, fill = Freq)) +
  # Dibujar cada celda de la matriz
  geom_tile(color = "white") +
  
  # Mostrar la frecuencia correspondiente en cada celda
  geom_text(aes(label = Freq), size = 6) +
  
  # Asignar un gradiente de color según la frecuencia
  scale_fill_gradient(low = "white", high = "steelblue") +
  
  # Invertir el eje Y para que la disposición coincida con la impresión de la matriz
  # de confusión generada por table() (la categoría 0 aparece en la fila superior)
  scale_y_discrete(limits = rev(levels(df$Real))) +
  
  # Añadir títulos a los ejes y a la leyenda
  labs(x = "Predicción",
       y = "Valor real",
       fill = "Frecuencia") +
  
  # Aplicar un tema minimalista al gráfico
  theme_minimal()


En un problema binario, la matriz se compone de cuatro cantidades fundamentales: verdaderos negativos (TN), verdaderos positivos (TP), falsos positivos (FP) y falsos negativos (FN).

Los cuatro casos posibles:

  • Verdaderos negativos (TN): el modelo predijo 0 y realmente era 0.
  • Verdaderos positivos (TP): el modelo predijo 1 y realmente era 1.
  • Falsos positivos (FP): el modelo predijo 1, pero realmente era 0 (una falsa alarma).
  • Falsos negativos (FN): el modelo predijo 0, pero realmente era 1 (se le escapó un caso).

A partir de esta tabla calculamos las métricas más usadas para evaluar un clasificador:

TP <- matriz_confusion["1", "1"]
TN <- matriz_confusion["0", "0"]
FP <- matriz_confusion["0", "1"]
FN <- matriz_confusion["1", "0"]

accuracy  <- (TP + TN) / sum(matriz_confusion)
precision <- TP / (TP + FP)
recall    <- TP / (TP + FN)
f1        <- 2 * precision * recall / (precision + recall)

cat("Accuracy:", round(accuracy, 3), "\n")
## Accuracy: 0.82
cat("Precision:", round(precision, 3), "\n")
## Precision: 0.714
cat("Recall:", round(recall, 3), "\n")
## Recall: 0.833
cat("F1 score:", round(f1, 3), "\n")
## F1 score: 0.769


¿Qué significa cada una?

  • Accuracy (exactitud): proporción de predicciones correctas sobre el total. Con nuestros datos, el modelo acierta en aproximadamente el 82% de los casos.
  • Precision (precisión): de todos los clientes a los que el modelo les dijo “sí van a comprar”, ¿qué porcentaje realmente compró? En nuestro caso, alrededor del 71%. Una precisión alta significa pocas falsas alarmas.
  • Recall (sensibilidad): de todos los clientes que realmente compraron, ¿a qué porcentaje detectó el modelo? Aproximadamente 83%. Un recall alto significa que se nos escapan pocos compradores reales.
  • F1 score: un balance entre precision y recall en un solo número (77% en nuestro caso). Es útil cuando nos importan ambas cosas por igual y no solo el accuracy.

Nota: el accuracy mide la proporción total de predicciones correctas, pero no siempre refleja la calidad real del modelo. Cuando una de las clases es mucho más frecuente que la otra, esta métrica puede resultar engañosa. Por ejemplo, si el 95 % de los clientes no compra un producto, un modelo que siempre prediga “no compra” alcanzaría un 95 % de accuracy, aunque sería incapaz de identificar a los clientes que sí compran. Por ello, es recomendable evaluar el modelo considerando en conjunto el accuracy, la precisión (precision), la sensibilidad (recall) y la especificidad (specificity), ya que cada una aporta información sobre un aspecto distinto de su desempeño.

Conclusión

Con ingreso mensual como único predictor, el modelo logra un desempeño razonable (82% de accuracy, buen balance entre precision y recall). No es perfecto ya que hay 12 falsos positivos y 6 falsos negativos sobre 100 clientes de prueba, pero es un punto de partida sólido, y coherente con lo que vimos en la exploración: el ingreso sí separa razonablemente bien a quienes compran de quienes no.

Ya vimos cómo abordar un problema de clasificación con dos categorías usando regresión logística. En la próxima sección veremos el análisis discriminante, otro enfoque para el mismo tipo de problema.

Análisis discriminante lineal (LDA)

El análisis discriminante lineal (LDA) también resuelve problemas de clasificación, pero partiendo de una idea distinta: en vez de modelar directamente la probabilidad de cada clase, LDA describe cómo se distribuyen los predictores dentro de cada clase, y usa eso para decidir a qué clase asignar una observación nueva.

Dos ventajas prácticas de LDA frente a la regresión logística, sin entrar en la demostración matemática:

El problema y los datos

Un equipo de biólogos quiere identificar a qué especie (a o b) pertenece un insecto, a partir de tres medidas: longitud de la pata, diámetro del abdomen y diámetro del órgano sexual. Se midieron 10 individuos de cada especie.

library(readxl)
library(ggplot2)
library(dplyr)
library(MASS) # ya viene instalada con R base, no hace falta install.packages()

insectos <- read_excel("Insectos.xlsx")

head(insectos)
## # A tibble: 6 × 4
##   Especie  Pata Abdomen Organo_Sexual
##   <chr>   <dbl>   <dbl>         <dbl>
## 1 a         191     131            53
## 2 a         185     134            50
## 3 a         200     137            52
## 4 a         173     127            50
## 5 a         171     128            49
## 6 a         160     118            47


insectos$Especie <- as.factor(insectos$Especie)

str(insectos)
## tibble [20 × 4] (S3: tbl_df/tbl/data.frame)
##  $ Especie      : Factor w/ 2 levels "a","b": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Pata         : num [1:20] 191 185 200 173 171 160 188 186 174 163 ...
##  $ Abdomen      : num [1:20] 131 134 137 127 128 118 134 129 131 115 ...
##  $ Organo_Sexual: num [1:20] 53 50 52 50 49 47 54 51 52 47 ...


Exploración rápida

Antes de modelar, veamos si las tres variables realmente diferencian a las especies.

ggplot(insectos, aes(x = Especie, y = Pata, fill = Especie)) +
  geom_boxplot() +
  labs(title = "Longitud de la pata según especie", x = "Especie", y = "Longitud de la pata") +
  theme(legend.position = "none")


Se ve una diferencia notoria: la especie b tiende a tener patas más largas, aunque hay algo de solapamiento entre 184 y 200. Podemos revisar las otras dos variables de la misma forma:

ggplot(insectos, aes(x = Especie, y = Abdomen, fill = Especie)) +
  geom_boxplot() +
  labs(title = "Diámetro del abdomen según especie", x = "Especie", y = "Abdomen") +
  theme(legend.position = "none")

ggplot(insectos, aes(x = Especie, y = Organo_Sexual, fill = Especie)) +
  geom_boxplot() +
  labs(title = "Diámetro del órgano sexual según especie", x = "Especie", y = "Órgano sexual") +
  theme(legend.position = "none")


Abdomen y Organo_Sexual muestran más solapamiento entre especies que Pata. Ninguna variable por sí sola separa perfectamente a las dos especies por lo que tiene sentido combinar las tres en un solo modelo.

LDA asume que, dentro de cada especie, los predictores se distribuyen aproximadamente como una campana (normal) y que la dispersión es parecida entre especies. Verificar esto formalmente requiere tests como Shapiro-Wilk o Box’s M, que no vamos a cubrir en el curso. Con un dataset tan chico (20 observaciones), esos tests tampoco son muy confiables. Nos alcanza con la revisión visual de arriba para seguir adelante.

Ajustando el modelo con lda()

La función lda() viene en el paquete MASS, que se instala junto con R (no hace falta install.packages(), solo library(MASS)).

modelo_lda <- lda(Especie ~ Pata + Abdomen + Organo_Sexual, data = insectos)

modelo_lda
## Call:
## lda(Especie ~ Pata + Abdomen + Organo_Sexual, data = insectos)
## 
## Prior probabilities of groups:
##   a   b 
## 0.5 0.5 
## 
## Group means:
##    Pata Abdomen Organo_Sexual
## a 179.1   128.4          50.5
## b 208.2   122.8          48.9
## 
## Coefficients of linear discriminants:
##                       LD1
## Pata           0.13225339
## Abdomen       -0.07941509
## Organo_Sexual -0.52655608


De esta salida nos interesan dos cosas:

  • Group means: el promedio de cada variable dentro de cada especie. Confirma lo que vimos en los boxplots. Por ejemplo, la especie b tiene en promedio patas más largas (208.2 vs. 179.1).
  • Coefficients of linear discriminants (LD1): son los “pesos” que LDA le da a cada variable para construir una única combinación (la función discriminante) que mejor separa las especies. El signo indica la dirección: Organo_Sexual tiene el coeficiente de mayor magnitud, lo que sugiere que es la variable que más aporta a separar las especies dentro de esta combinación — aunque, como vimos en los boxplots, no es la que más se distingue mirándola sola. Esto pasa gracias a que LDA combina las tres variables a la vez, no cada una por separado.

Prediciendo una observación nueva

Supongamos que llega un nuevo insecto con pata = 194, abdomen = 124 y órgano sexual = 49. ¿A qué especie lo asignamos?

nuevo_insecto <- data.frame(Pata = 194, Abdomen = 124, Organo_Sexual = 49)

predict(modelo_lda, newdata = nuevo_insecto)
## $class
## [1] b
## Levels: a b
## 
## $posterior
##            a         b
## 1 0.05823333 0.9417667
## 
## $x
##         LD1
## 1 0.5419421


El resultado devuelve tres partes:

  • class: la especie asignada (la de mayor probabilidad).
  • posterior: la probabilidad estimada de pertenecer a cada especie.
  • x: el valor de la función discriminante para esta observación.

Para este insecto, el modelo predice la especie b con una probabilidad de aproximadamente 94%.

Evaluando el modelo: train / test

Como en las secciones anteriores, para evaluar el modelo separamos entrenamiento y prueba.

set.seed(123)

n <- nrow(insectos)
indices_entrenamiento <- sample(1:n, size = 0.8 * n)

train <- insectos[indices_entrenamiento, ]
test <- insectos[-indices_entrenamiento, ]

nrow(train)
## [1] 16
nrow(test)
## [1] 4


Con solo 20 observaciones en total, un split 80/20 deja apenas 4 observaciones de prueba. Es un dataset de juguete pensado para enseñar el método, no para sacar conclusiones generales pero suficiente para entender el flujo de trabajo.

modelo_train <- lda(Especie ~ Pata + Abdomen + Organo_Sexual, data = train)

pred_test <- predict(modelo_train, newdata = test)

matriz_confusion <- table(Real = test$Especie, Predicho = pred_test$class)
matriz_confusion
##     Predicho
## Real a b
##    a 2 0
##    b 1 1


Calculamos las mismas métricas que en regresión logística:

TP <- matriz_confusion["b", "b"]
TN <- matriz_confusion["a", "a"]
FP <- matriz_confusion["a", "b"]
FN <- matriz_confusion["b", "a"]

accuracy  <- (TP + TN) / sum(matriz_confusion)
precision <- TP / (TP + FP)
recall    <- TP / (TP + FN)

# Media armónica entre precision y recall
f1 <- 2 * (precision * recall) / (precision + recall)

cat("Accuracy :", round(accuracy, 3), "\n")
## Accuracy : 0.75
cat("Precision:", round(precision, 3), "\n")
## Precision: 1
cat("Recall   :", round(recall, 3), "\n")
## Recall   : 0.5
cat("F1-score :", round(f1, 3), "\n")
## F1-score : 0.667


Con la semilla 123, el modelo clasifica correctamente 3 de las 4 observaciones de prueba (75% de accuracy): se equivoca en un insecto de la especie b que quedó clasificado como a.

Sin embargo, este resultado debe interpretarse con cautela. Al contar con solo 4 observaciones de prueba, un único error modifica el accuracy en 25 puntos porcentuales. Es la misma advertencia que vimos en regresión lineal con datasets pequeños: las métricas de evaluación pueden variar considerablemente por el reducido tamaño de la muestra. En consecuencia, con un conjunto de datos tan pequeño no conviene extraer conclusiones firmes sobre el desempeño del modelo. En este ejemplo, el objetivo principal es comprender el flujo de trabajo para aplicar un análisis discriminante lineal (LDA), más que evaluar su capacidad predictiva.

Tanto la regresión logística como el LDA resuelven problemas de clasificación partiendo de predictores que describen cada observación. En la próxima sección cambiamos de objetivo para en vez de predecir una categoría, encontrar cómo reducir un conjunto grande de variables a unas pocas dimensiones que concentren la mayor parte de la información a través del Análisis de Componentes Principales (PCA).

Análisis de Componentes Principales (PCA)

Hasta ahora, todos los modelos que vimos en el curso (regresión lineal, regresión logística y LDA) comparten algo en común: en todos existía una variable respuesta Y que queríamos predecir a partir de un conjunto de predictores X. A este tipo de problemas se les llama aprendizaje supervisado.

Hoy vamos a estudiar un problema distinto. El Análisis de Componentes Principales (Principal Component Analysis, o PCA) no busca predecir nada. No hay una variable Y, solo tenemos un conjunto de variables X y queremos entender mejor su estructura.

Por eso decimos que PCA es un método de aprendizaje no supervisado: no hay una “respuesta correcta” contra la cual comparar los resultados, el objetivo es puramente descubrir patrones dentro de los datos mismos.

Concretamente, PCA sirve para responder preguntas como:

Esta última pregunta va a ser justamente con la que hagamos cierre de la clase de hoy.

La idea intuitiva: ¿qué hace PCA?

Un ejemplo de juguete en 2 dimensiones

Antes de usar un dataset real, construyamos un ejemplo simple para entender la idea central de PCA: dos variables inventadas y correlacionadas entre sí.

set.seed(123)

x <- rnorm(60, mean = 50, sd = 10)
y <- 0.8 * x + rnorm(60, mean = 0, sd = 5)

juguete <- data.frame(x, y)

plot(juguete$x, juguete$y,
     pch = 19, col = "steelblue",
     xlab = "Variable 1", ylab = "Variable 2",
     main = "Dos variables correlacionadas")


La nube de puntos no es circular, tiene una forma alargada, como un “cigarro” inclinado. Eso es justamente lo que significa que dos variables estén correlacionadas: si conocemos el valor de x, podemos tener una idea de cuánto vale y.

¿Qué hace PCA con esta nube de puntos? Busca un nuevo eje, una nueva dirección, que capture la mayor cantidad de variación posible en los datos. Después busca un segundo eje, perpendicular al primero, que capture la variación que sobró.

# PCA sobre los datos estandarizados (media = 0, desvío estándar = 1)
pca_juguete <- prcomp(juguete, scale. = TRUE)

# Obtener explícitamente los datos estandarizados
juguete_std <- scale(juguete)

# Como los datos están estandarizados, el centro es aproximadamente (0, 0)
centro <- colMeans(juguete_std)

# Graficar los datos estandarizados
# asp = 1 asegura que la misma distancia en X e Y tenga el mismo tamaño
# en la pantalla, de modo que las direcciones ortogonales se vean perpendiculares.
plot(juguete_std[, 1], juguete_std[, 2],
     asp = 1,
     pch = 19, col = "steelblue",
     xlab = "Variable 1 (estandarizada)",
     ylab = "Variable 2 (estandarizada)",
     main = "Direcciones encontradas por PCA")

# Dirección del primer componente principal (máxima varianza)
arrows(centro[1], centro[2],
       centro[1] + pca_juguete$rotation[1, 1] * 3,
       centro[2] + pca_juguete$rotation[2, 1] * 3,
       col = "red", lwd = 2, length = 0.15)

# Dirección del segundo componente principal
# Es perpendicular al primero en el espacio estandarizado.
arrows(centro[1], centro[2],
       centro[1] + pca_juguete$rotation[1, 2] * 3,
       centro[2] + pca_juguete$rotation[2, 2] * 3,
       col = "darkgreen", lwd = 2, length = 0.15)

legend("topleft",
       legend = c("Componente 1", "Componente 2"),
       col = c("red", "darkgreen"),
       lwd = 2,
       bty = "n")


#--------------------------------------------------------------------------
# 1. Calcular el Análisis de Componentes Principales (PCA)
#--------------------------------------------------------------------------

# scale = TRUE indica que antes de calcular el PCA cada variable se
# estandariza:
#   - se resta su media
#   - se divide por su desvíación estándar
#
# De esta forma todas las variables quedan con:
#   media = 0
#   desvío estándar = 1
#
# Esto evita que una variable "domine" el análisis simplemente porque
# tiene valores que en magnitud numérica son más grandes.
pca_juguete <- prcomp(juguete, scale. = TRUE)


#--------------------------------------------------------------------------
# 2. Obtener los datos estandarizados
#--------------------------------------------------------------------------

# scale() devuelve exactamente los datos utilizados por prcomp().
#
# Es importante graficar estos datos y no los originales, ya que las
# direcciones (componentes principales) fueron calculadas sobre este
# espacio estandarizado.
juguete_std <- scale(juguete)


#--------------------------------------------------------------------------
# 3. Calcular el centro del gráfico
#--------------------------------------------------------------------------

# Después de estandarizar, la media de cada variable es aproximadamente 0.
#
# Es decir, el "centro" de la nube de puntos queda ubicado en (0,0).
#
# Desde este punto se dibujarán las flechas que representan las componentes
# principales.
centro <- colMeans(juguete_std)


#--------------------------------------------------------------------------
# 4. Dibujar la nube de puntos
#--------------------------------------------------------------------------

plot(juguete_std[,1], juguete_std[,2],

     # asp = 1 obliga a que una unidad en X tenga exactamente el mismo
     # tamaño que una unidad en Y.
     #
     # Sin esta opción, un ángulo recto podría verse inclinado debido a la
     # forma del gráfico.
     asp = 1,

     # pch = 19 dibuja puntos sólidos.
     pch = 19,

     # Color de los puntos.
     col = "steelblue",

     # Etiquetas de los ejes.
     xlab = "Variable 1 (estandarizada)",
     ylab = "Variable 2 (estandarizada)",

     # Título del gráfico.
     main = "Direcciones encontradas por PCA"
)


#--------------------------------------------------------------------------
# 5. Dibujar la primera componente principal
#--------------------------------------------------------------------------

# La matriz rotation contiene la dirección de cada componente principal.
#
# La primera columna corresponde a la primera componente (PC1), que es la
# dirección donde los datos presentan la mayor variabilidad.
#
# arrows(x0, y0, x1, y1) dibuja una flecha desde el punto inicial
# (x0,y0) hasta el punto final (x1,y1).

arrows(
  centro[1], centro[2],                   # Punto inicial (el centro)

  # Coordenada X del extremo de la flecha.
  # Multiplicamos por 3 únicamente para hacer la flecha más visible.
  centro[1] + pca_juguete$rotation[1,1] * 3,

  # Coordenada Y del extremo.
  centro[2] + pca_juguete$rotation[2,1] * 3,

  # Apariencia de la flecha.
  col = "red",
  lwd = 2,
  length = 0.15
)


#--------------------------------------------------------------------------
# 6. Dibujar la segunda componente principal
#--------------------------------------------------------------------------

# La segunda columna de rotation contiene la segunda componente (PC2).
#
# Esta componente siempre es perpendicular a la primera y explica la mayor
# cantidad posible de variabilidad restante.
#
# Juntas forman un nuevo sistema de ejes para describir los datos.

arrows(
  centro[1], centro[2],

  centro[1] + pca_juguete$rotation[1,2] * 3,
  centro[2] + pca_juguete$rotation[2,2] * 3,

  col = "darkgreen",
  lwd = 2,
  length = 0.15
)


#--------------------------------------------------------------------------
# 7. Agregar una leyenda
#--------------------------------------------------------------------------

# La leyenda identifica qué flecha representa cada componente principal.
legend(
  "topleft",
  legend = c("Componente 1", "Componente 2"),
  col = c("red", "darkgreen"),
  lwd = 2,
  bty = "n"   # Sin borde alrededor de la leyenda
)

Nota: el factor por 3 solo hace que las flechas sean más largas para facilitar su visualización. La información relevante es su dirección, no su longitud.

La flecha roja (Componente Principal 1, o PC1) apunta justo a lo largo del “cigarro”: es la dirección donde los datos varían más. La flecha verde (Componente Principal 2, o PC2) es perpendicular a la primera, y apunta hacia donde queda la variación restante, que en este ejemplo es poca, porque la nube es angosta en esa dirección.

Esa es la idea detrás de PCA:

  • Cada componente principal es una combinación lineal de las variables originales.
  • El primer componente es la dirección de máxima varianza posible.
  • El segundo componente es la dirección de máxima varianza posible que sea perpendicular (no correlacionada) con el primero. Y así sucesivamente para el tercero, cuarto, etc.
  • Cuando tenemos p variables originales, PCA siempre nos entrega p componentes principales, pero la esperanza es que unos pocos ya capturen casi toda la información, y podamos ignorar el resto.

No vamos a profundizar en el álgebra lineal detrás de cómo se calculan estas direcciones (eso involucra descomposición en autovalores/autovectores). Para este curso nos basta con entender qué representa cada pieza y cómo interpretarla, ya que R se encarga de todo el cálculo por nosotros.

Trabajando con datos reales: USArrests

USArrests es un dataset clásico que ya viene incluido en R (no hace falta descargar nada). Contiene, para cada uno de los 50 estados de Estados Unidos, la tasa de arrestos por cada 100,000 habitantes en tres delitos (asesinato, asalto y violación), además del porcentaje de población que vive en zonas urbanas.

Explorando el dataset

datos <- USArrests

str(datos)
## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...


head(datos)
##            Murder Assault UrbanPop Rape
## Alabama      13.2     236       58 21.2
## Alaska       10.0     263       48 44.5
## Arizona       8.1     294       80 31.0
## Arkansas      8.8     190       50 19.5
## California    9.0     276       91 40.6
## Colorado      7.9     204       78 38.7


Tenemos 50 observaciones (los estados quedan como nombres de fila, no como columna) y 4 variables numéricas: Murder, Assault, UrbanPop y Rape.

Como UrbanPop registra un porcentaje, un registro de 80 nos dice que el 80% de los habitantes de este estado vive en áreas urbanas, mientras que el 20% restante vive en áreas rurales. Esta característica demográfica del estado es considerada porque el grado de urbanización podría estar relacionado con los niveles de ciertos delitos.

El problema de la redundancia

Antes de aplicar PCA, veamos qué tan correlacionadas están estas 4 variables entre sí.

round(cor(datos), 2)
##          Murder Assault UrbanPop Rape
## Murder     1.00    0.80     0.07 0.56
## Assault    0.80    1.00     0.26 0.67
## UrbanPop   0.07    0.26     1.00 0.41
## Rape       0.56    0.67     0.41 1.00


Murder, Assault y Rape están bastante correlacionadas entre sí (todas rondan 0.6-0.8). Esto tiene sentido: un estado con más criminalidad en general probablemente tenga tasas más altas en varios tipos de delito a la vez. UrbanPop es la que menos se relaciona con las demás.

Podemos verlo también de forma visual con pairs(), que arma automáticamente un gráfico de dispersión para cada par de variables.

pairs(datos, pch = 19, col = "steelblue")


Esta correlación es justamente el síntoma que hace útil a PCA: si varias variables “cuentan la misma historia” (varían juntas), probablemente no necesitemos las 4 variables originales para explicar casi toda la información contenida en el dataset.

¿Por qué es importante escalar los datos?

sapply(datos, sd)
##    Murder   Assault  UrbanPop      Rape 
##  4.355510 83.337661 14.474763  9.366385


Assault tiene una desviación estándar (~85) enormemente más grande que las demás variables. PCA busca direcciones de máxima varianza, así que si no hacemos nada al respecto, Assault dominaría completamente el resultado simplemente por estar medida en una escala más grande, y no porque sea “más importante”.

La solución es estandarizar cada variable (restar su media y dividir entre su desviación estándar) antes de aplicar PCA, para que todas partan de la misma escala. En R, esto se controla con el argumento scale. de prcomp().

Regla práctica: salvo que todas tus variables ya estén en las mismas unidades y en escalas comparables, casi siempre conviene usar scale. = TRUE.

Ajustando PCA en R: prcomp()

pca_arrestos <- prcomp(datos, scale. = TRUE)

pca_arrestos
## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995 -0.4181809  0.3412327  0.64922780
## Assault  -0.5831836 -0.1879856  0.2681484 -0.74340748
## UrbanPop -0.2781909  0.8728062  0.3780158  0.13387773
## Rape     -0.5434321  0.1673186 -0.8177779  0.08902432


El resultado tiene dos piezas centrales que vamos a revisar por separado: rotation y x.

Los “loadings” o cargas: rotation

pca_arrestos$rotation
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995 -0.4181809  0.3412327  0.64922780
## Assault  -0.5831836 -0.1879856  0.2681484 -0.74340748
## UrbanPop -0.2781909  0.8728062  0.3780158  0.13387773
## Rape     -0.5434321  0.1673186 -0.8177779  0.08902432


Esta matriz nos dice cómo se construye cada componente principal a partir de las variables originales. Por ejemplo, la columna PC1 son los “pesos” con los que se combinan Murder, Assault, UrbanPop y Rape para formar el primer componente.

Podemos interpretar el signo y la magnitud de estos pesos:

  • En PC1, los pesos de Murder, Assault y Rape son similares en magnitud y del mismo signo, mientras que UrbanPop tiene un peso mucho menor. Esto sugiere que PC1 resume principalmente el nivel general de criminalidad violenta del estado.
  • En PC2, el patrón se invierte: UrbanPop domina. Este componente separa a los estados principalmente por qué tan urbanizados son, casi sin importar su nivel de criminalidad.

Esta es la gran ventaja interpretativa de PCA: no solo reduce dimensiones, también puede revelarnos qué “historias” subyacentes explican la variación en los datos.

Los “scores” o nuevas coordenadas: x

head(pca_arrestos$x)
##                   PC1        PC2         PC3          PC4
## Alabama    -0.9756604 -1.1220012  0.43980366  0.154696581
## Alaska     -1.9305379 -1.0624269 -2.01950027 -0.434175454
## Arizona    -1.7454429  0.7384595 -0.05423025 -0.826264240
## Arkansas    0.1399989 -1.1085423 -0.11342217 -0.180973554
## California -2.4986128  1.5274267 -0.59254100 -0.338559240
## Colorado   -1.4993407  0.9776297 -1.08400162  0.001450164


Mientras que rotation nos dice cómo se construyen los componentes, x contiene el resultado de aplicar esa combinación a cada observación: es decir, la posición de cada estado en el nuevo sistema de coordenadas formado por los componentes principales.

Por ejemplo, un estado con un valor muy alto en PC1 es un estado con criminalidad violenta muy por encima del promedio.

Visualizando con biplot()

R nos permite ver rotation y x al mismo tiempo con una sola función base, biplot().

biplot(pca_arrestos, scale = 0, cex = 0.6)


Este gráfico combina dos cosas:

  • Los nombres de los estados (en negro) están ubicados según sus valores en PC1 (eje horizontal) y PC2 (eje vertical). Estados cercanos entre sí tienen un perfil de criminalidad/urbanización parecido.
  • Las flechas rojas representan las variables originales. La dirección de la flecha indica hacia dónde “empuja” esa variable, y su longitud indica qué tan fuerte es su contribución.

Por ejemplo, vemos que Murder, Assault y Rape apuntan casi en la misma dirección (hacia la derecha), confirmando que están correlacionadas y todas “empujan” en el sentido de PC1. UrbanPop apunta hacia arriba, en la dirección de PC2, separada de las demás.

Estados como California o Nevada, ubicados hacia la izquierda, tienen niveles de criminalidad violenta relativamente altos. Estados como Dakota del Norte, hacia la derecha, tienen niveles bajos.

¿Cuánta información conserva cada componente?

Con 4 variables originales, obtuvimos 4 componentes. Pero, ¿realmente necesitamos los 4? Para responder eso necesitamos cuantificar cuánta información (varianza) aporta cada uno.

summary(pca_arrestos)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion  0.6201 0.8675 0.95664 1.00000


Nos interesan dos filas de esta tabla:

  • Proportion of Variance: qué porcentaje de la variabilidad total del dataset es capturado por cada componente, por separado.
  • Cumulative Proportion: la suma acumulada de la fila anterior, es decir, cuánta varianza capturamos si usamos ese componente y todos los anteriores juntos.

Podemos extraer estos valores directamente para trabajar con ellos:

varianza_explicada <- summary(pca_arrestos)$importance[2, ] # fila "Proportion of Variance"
varianza_acumulada <- summary(pca_arrestos)$importance[3, ] # fila "Cumulative Proportion"

varianza_explicada
##     PC1     PC2     PC3     PC4 
## 0.62006 0.24744 0.08914 0.04336
varianza_acumulada
##     PC1     PC2     PC3     PC4 
## 0.62006 0.86750 0.95664 1.00000


El primer componente por sí solo explica alrededor del 62% de toda la variabilidad del dataset, y entre PC1 y PC2 juntos ya cubrimos cerca del 87%. Es decir, podríamos resumir un dataset de 4 variables en solamente 2 componentes y quedarnos con la gran mayoría de la información original.

Visualizando con un screeplot

Un screeplot (o “gráfico de codo”) es simplemente la varianza explicada de cada componente, graficada en orden. prcomp trae un método de plot() que lo hace directamente:

plot(pca_arrestos, type = "l", main = "Varianza explicada por componente")


También podemos armar la versión de varianza acumulada nosotros mismos, que suele ser más útil para decidir cuántos componentes conservar:

barplot(varianza_acumulada,
        names.arg = paste0("PC", 1:4),
        col = "steelblue",
        ylim = c(0, 1),
        ylab = "Varianza acumulada",
        main = "Varianza acumulada por componente")
abline(h = 0.9, col = "red", lty = 2) # referencia visual del 90%


¿Cuántos componentes principales elegir?

Esta es la pregunta más frecuente al usar PCA, y la respuesta honesta es: no existe una regla única ni exacta. Es una decisión que combina el gráfico anterior con el objetivo práctico que tengamos. Algunos criterios comunes que se usan en la práctica:

  • Regla del codo (elbow rule): en el screeplot, buscamos el punto donde la curva deja de caer pronunciadamente y se “aplana”. Los componentes antes del codo suelen concentrar la información importante; los de después, capturan mayormente ruido.

  • Varianza acumulada objetivo: fijar de antemano un porcentaje deseado (por ejemplo, 80% o 90%) y quedarnos con los primeros componentes que alcancen ese umbral. El umbral en sí también es una elección subjetiva: depende de qué tanto estemos dispuestos a sacrificar información por simplicidad.

  • Regla de Kaiser: conservar únicamente los componentes cuya varianza (el cuadrado de pca$sdev) sea mayor a 1 cuando trabajamos con datos escalados. La lógica es que, al estandarizar, cada variable original aporta exactamente 1 unidad de varianza, así que un componente “vale la pena” solo si logra capturar más varianza que una variable individual cualquiera.

pca_arrestos$sdev^2 # varianza de cada componente
## [1] 2.4802416 0.9897652 0.3565632 0.1734301


  • El objetivo del análisis: si el propósito es únicamente visualizar los datos en un gráfico, casi siempre nos quedamos con 2 (o a lo más 3) componentes, sin importar cuánta varianza acumulen exactamente. Si el propósito es usar los componentes como entrada de otro modelo (como haremos a continuación), conviene ser un poco más generosos y quedarnos con más componentes, para no perder información relevante para la predicción.

En resumen: el número de componentes no es un resultado que R nos entregue automáticamente como “correcto”, es una decisión que nosotros tomamos apoyándonos en estas herramientas.

Cerrando con un caso más exigente: cuando “más variables” no es mejor

Todo lo anterior lo vimos con un dataset chico y amigable (4 variables). Ahora vamos a construir un escenario deliberadamente más difícil, parecido a uno que podrían encontrarse en la práctica: un dataset con muchas variables, varias de ellas redundantes entre sí, y relativamente pocas observaciones.

Construyendo un dataset más complejo

Vamos a simular un dataset con 50 variables predictoras y 120 observaciones. La clave (que en la vida real nunca conoceríamos de antemano) es que, por diseño, toda esa información en realidad proviene de únicamente 3 factores subyacentes: las 50 variables no son más que versiones ruidosas y redundantes de esos 3 factores, medidos una y otra vez de formas ligeramente distintas.

set.seed(123)

n <- 120  # observaciones
p <- 50   # variables predictoras observadas
k <- 3    # factores "verdaderos" que generan los datos (en la práctica, desconocido)

# 3 factores subyacentes, no observables directamente
factores <- matrix(rnorm(n * k), nrow = n, ncol = k)

# Cada una de las 50 variables observadas es una combinación ruidosa de esos 3 factores
cargas <- matrix(runif(p * k, -1, 1), nrow = p, ncol = k)
ruido_x <- matrix(rnorm(n * p, sd = 0.5), nrow = n, ncol = p)

X <- factores %*% t(cargas) + ruido_x
colnames(X) <- paste0("var", 1:p)

# La variable respuesta depende únicamente de los 3 factores reales, más ruido propio
y <- 3 * factores[, 1] - 2 * factores[, 2] + 1.5 * factores[, 3] + rnorm(n, sd = 3)

datos_complejos <- data.frame(y = y, X)

dim(datos_complejos)
## [1] 120  51


head(datos_complejos[, 1:6])
##            y       var1       var2       var3      var4        var5
## 1 -3.1454841 -0.6473229 -0.6023780 -1.3778788 -0.709707 -0.04371788
## 2  0.4523749 -0.5840313 -0.5297276 -1.1109577 -1.052958  1.13369822
## 3  6.7188419  0.4453915 -0.4561263  0.2355781  2.104250  0.46630690
## 4 -2.7552533 -1.7662910  0.2010061 -0.7228665 -1.106391 -0.14110012
## 5 -0.7905201  2.0632777  0.9638096  0.3268852  1.378077 -0.63411124
## 6  8.8722264 -0.7238704 -0.7958800  1.1954791  2.921448  1.57328742


Hasta ahora usamos PCA para explorar y resumir un conjunto de variables sin intentar predecir nada. A partir de ahora el objetivo será construir un modelo de regresión para predecir y.

Nota: PCA continúa siendo una técnica no supervisada: se aplica exclusivamente sobre las variables predictoras (X) para obtener un conjunto reducido de componentes, que luego se utilizarán como predictores del modelo.

Revisemos la correlación entre estas variables.

round(cor(datos_complejos[, 2:7]), 2)
##       var1  var2  var3  var4  var5  var6
## var1  1.00  0.32  0.38  0.49 -0.42 -0.60
## var2  0.32  1.00  0.42  0.20 -0.62  0.23
## var3  0.38  0.42  1.00  0.56 -0.21 -0.19
## var4  0.49  0.20  0.56  1.00  0.10 -0.49
## var5 -0.42 -0.62 -0.21  0.10  1.00 -0.08
## var6 -0.60  0.23 -0.19 -0.49 -0.08  1.00


Muchas correlaciones altas entre variables que, en teoría, “deberían” ser predictores independientes.

Regresión lineal múltiple con todos los predictores

Apliquemos lo que ya sabemos: ajustemos un modelo de regresión lineal múltiple usando las 50 variables como predictores. La sintaxis y ~ . le dice a R “usa todas las demás columnas del data frame como predictores”.

modelo_completo <- lm(y ~ ., data = datos_complejos)

summary(modelo_completo)$r.squared
## [1] 0.7720167
summary(modelo_completo)$adj.r.squared
## [1] 0.6068113


Aquí es donde vale la pena detenerse. Es muy probable que encuentren una brecha bastante grande entre el R² normal y el R² ajustado: el R² puede verse razonablemente alto, pero el R² ajustado cae de forma notoria.

Como usamos 50 predictores para explicar y, pero muchos de ellos son versiones redundantes entre sí de la misma información, el R² ajustado está penalizando fuertemente el hecho de que estemos “gastando” 50 grados de libertad para, en el fondo, aprovechar mucha menos información genuina que esa cantidad de variables sugeriría.

Con tan pocas observaciones (120) para tantos predictores (50), el modelo tiene mucho margen para ajustarse a particularidades del propio dataset de entrenamiento, más que a una relación real y generalizable.

La sospecha de que el modelo se está ajustando a particularidades de estos datos puntuales puede verificarse dividiendo los datos en entrenamiento y prueba.

Evaluando la generalización: entrenamiento vs. prueba

set.seed(123)

n_total <- nrow(datos_complejos)
indices_entrenamiento <- sample(1:n_total, size = 0.75 * n_total)

train <- datos_complejos[indices_entrenamiento, ]
test <- datos_complejos[-indices_entrenamiento, ]

nrow(train)
## [1] 90
nrow(test)
## [1] 30


modelo_completo_train <- lm(y ~ ., data = train)

# Error sobre los datos de entrenamiento
pred_train <- predict(modelo_completo_train, newdata = train)
rmse_train_completo <- sqrt(mean((train$y - pred_train)^2))

# Error sobre datos que el modelo nunca vio
pred_test <- predict(modelo_completo_train, newdata = test)
rmse_test_completo <- sqrt(mean((test$y - pred_test)^2))

rmse_train_completo
## [1] 1.938255
rmse_test_completo
## [1] 3.748312


Con 50 predictores y solo 90 observaciones de entrenamiento (75% de 120), el modelo tiene muchísima flexibilidad para ajustarse casi a la perfección a los datos de entrenamiento. El problema es que buena parte de ese ajuste es puro sobreajuste: el error crece de forma marcada al evaluar sobre test, datos que el modelo nunca vio. Esta es la señal más contundente de que el modelo completo no está generalizando bien, más contundente incluso que la brecha entre R² y R² ajustado que vimos antes.

Aplicando PCA para reducir dimensionalidad

En lugar de usar las 50 variables originales (correlacionadas y ruidosas) como predictores, probemos resumirlas primero en unos pocos componentes principales, y usemos esos componentes como los nuevos predictores.

Importante: el PCA lo ajustamos únicamente con train, igual que hicimos con el modelo de regresión. Si usáramos también test para calcular los componentes, estaríamos “espiando” información que se supone el modelo no debería ver todavía.

# Excluimos la columna "y" porque PCA se aplica solo sobre los predictores
pca_train <- prcomp(train[, -1], scale. = TRUE)

summary(pca_train)$importance[, 1:6]
##                             PC1      PC2      PC3      PC4      PC5       PC6
## Standard deviation     3.736805 3.529133 3.333253 1.025307 0.919325 0.8690302
## Proportion of Variance 0.279270 0.249100 0.222210 0.021030 0.016900 0.0151000
## Cumulative Proportion  0.279270 0.528370 0.750580 0.771610 0.788510 0.8036100


plot(pca_train, type = "l", main = "Varianza explicada por componente")


Debería verse un “codo” bastante marcado: los primeros componentes concentran una parte grande de la varianza, y después de un punto, cada componente adicional aporta muy poco. Esto es exactamente lo que esperábamos, porque construimos el dataset a partir de solamente 3 factores reales.

Con base en ese codo, nos quedamos con los primeros 3 componentes.

varianza_acumulada_compleja <- summary(pca_train)$importance[3, ]
varianza_acumulada_compleja[1:5]
##     PC1     PC2     PC3     PC4     PC5 
## 0.27927 0.52837 0.75058 0.77161 0.78851


Con solamente 3 de los 50 componentes posibles, ya estamos reteniendo una fracción sustancial de toda la variabilidad original del dataset. Pasamos de 50 variables a 3, sin haber mirado ni una sola vez la variable y en el proceso.

Regresión con los componentes principales

Armamos un nuevo data frame de entrenamiento, reemplazando las 50 variables originales por los 3 componentes principales, y ajustamos la regresión sobre ese data frame reducido.

scores_train <- as.data.frame(pca_train$x[, 1:3])
scores_train$y <- train$y

modelo_pca <- lm(y ~ ., data = scores_train)

summary(modelo_pca)$r.squared
## [1] 0.663538
summary(modelo_pca)$adj.r.squared
## [1] 0.651801


Ahora la brecha entre R² y R² ajustado es mucho más pequeña. Ahora solo estamos “pagando” el costo de 3 predictores en lugar de 50, así que la penalización por grados de libertad es mínima.

Para evaluar qué tan bien generaliza este modelo, necesitamos proyectar también los datos de test sobre los mismos 3 componentes que aprendimos con train. La función predict() funciona también sobre objetos de tipo prcomp, y logra aplicar la misma transformación (misma media, misma desviación estándar, mismos pesos de rotation) a datos nuevos.

scores_test <- as.data.frame(predict(pca_train, newdata = test[, -1])[, 1:3])
scores_test$y <- test$y

pred_train_pca <- predict(modelo_pca, newdata = scores_train)
rmse_train_pca <- sqrt(mean((scores_train$y - pred_train_pca)^2))

pred_test_pca <- predict(modelo_pca, newdata = scores_test)
rmse_test_pca <- sqrt(mean((scores_test$y - pred_test_pca)^2))

rmse_train_pca
## [1] 2.682392
rmse_test_pca
## [1] 2.731542


Resumen comparativo de ambos enfoques

data.frame(
  modelo = c("Regresión con las 50 variables originales", "Regresión con 3 componentes principales"),
  R2_ajustado = c(summary(modelo_completo_train)$adj.r.squared, summary(modelo_pca)$adj.r.squared),
  RMSE_train = c(rmse_train_completo, rmse_train_pca),
  RMSE_test = c(rmse_test_completo, rmse_test_pca)
)
##                                      modelo R2_ajustado RMSE_train RMSE_test
## 1 Regresión con las 50 variables originales   0.5990974   1.938255  3.748312
## 2   Regresión con 3 componentes principales   0.6518010   2.682392  2.731542


Esta tabla resume toda la lección:

  • El modelo con las 50 variables originales probablemente tenga un RMSE de entrenamiento muy bajo (se ajusta casi perfectamente a esos datos puntuales), pero un RMSE de prueba notoriamente más alto, lo que exhibe sobreajuste.
  • El modelo con solo 3 componentes principales tiene un RMSE de entrenamiento un poco más alto que el modelo completo (tiene menos flexibilidad, es “menos poderoso” para memorizar), pero un RMSE de prueba mucho más bajo y parecido a su propio RMSE de entrenamiento, lo que nos dice que el modelo generaliza razonablemente bien.

¿Por qué ocurre esto? Porque las 50 variables originales, al ser tan redundantes entre sí, no le estaban dando al modelo completo 50 fuentes de información independiente, sino apenas 3 “reales” repetidas de muchas formas ligeramente distintas, más un montón de ruido individual de cada variable. El modelo completo termina usando ese exceso de parámetros para ajustarse al ruido específico del conjunto de entrenamiento. PCA, al resumir la información redundante en unos pocos componentes, elimina gran parte de ese ruido antes de que llegue siquiera al modelo de regresión.

Además de la ventaja en generalización, el modelo con componentes principales tiene otras ventajas prácticas:

  • Es un modelo mucho más simple de reportar e interpretar: 3 coeficientes en lugar de 50.
  • Al ser los componentes principales no correlacionados entre sí por construcción, no hay problema de multicolinealidad entre los predictores del modelo final, algo que sí podía ser un problema real en el modelo con las 50 variables originales correlacionadas.

La contraparte, y por eso PCA no es una solución mágica ni gratuita, es que perdimos algo de interpretabilidad directa: ya no tenemos coeficientes para “var1”, “var2”, etc., sino para combinaciones abstractas de todas ellas. Ese es exactamente el tipo de trade-off (simplicidad y generalización vs. interpretabilidad directa) que hay que sopesar cada vez que decidamos aplicar PCA como paso previo a otro modelo (como una regresión).