Clase 5 · Regresión lineal y regularización

Analítica de Datos · Maestría en Ciencias del Comportamiento · Universidad de San Andrés

Primavera 2026 · 05/09/2026

Abrir en Colab

Esta es la versión en R de la notebook de la clase. Cubre exactamente los mismos pasos, los mismos datos y los mismos números que la de Python: si seguiste una, la otra no trae contenido nuevo.

La pregunta de hoy es ¿cuánto bienestar laboral podemos predecir a partir del salario?, y de ahí en adelante, ¿cuánto mejora si usamos más variables?

El modelo que mejor ajusta los datos que ya viste no es el que mejor predice los que vienen.

# Paso Qué hacemos
1 Armar la tabla unir nimbus_clima con nimbus_salario
2 Una recta a ojo mover β₀ y β₁ a mano y mirar los residuos
3 Mínimos cuadrados programar la fórmula, sin lm()
4 Medir el ajuste RSE y R², a mano y después con lm()
5 Varias variables ir agregando predictores y ver qué pasa
6 El número honesto partir en entrenamiento y testeo, y volver a medir

Ojo con los dos “bienestar”. El de hoy es bienestar_laboral, un índice de 0 a 100 de la encuesta de clima 2026. No es el bienestar diario del piloto de la fruta (escala 1 a 7).

1. Armar la tabla

clima tiene una fila por empleado. salario tiene una fila por empleado y por año, así que primero hay que quedarse con un año.

library(dplyr)
library(ggplot2)

# Todo sale del espejo público de la materia. Si cambia de lugar, se toca esta línea.
BASE <- "https://raw.githubusercontent.com/tomdamelio/analitica_de_datos_alumnos/main/data/toy-nimbus/"

clima <- read.csv(paste0(BASE, "nimbus_clima.csv"))
salario <- read.csv(paste0(BASE, "nimbus_salario.csv"))

cat("clima  ", dim(clima), "\n")
cat("salario", dim(salario), "\n")
head(clima, 3)

Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


clima   600 20 
salario 1800 4
A data.frame: 3 × 20
empleado_id bienestar_laboral horas_extra_semana apoyo_equipo reconocimiento autonomia dias_home_office bono_anual_pct reuniones_semana mensajes_chat_dia dias_vacaciones_tomados distancia_oficina_km cursos_completados meses_en_el_rol tickets_cerrados_mes emails_enviados_dia proyectos_activos dias_licencia_medica horas_capacitacion puntualidad_pct
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <int> <int> <int> <dbl> <int> <int> <int> <int> <int> <int> <dbl> <dbl>
1 1 66.2 3.2 7.7 5.1 7.4 0 11.31 9 43 10 8.9 1 32 15 11 2 4 7.4 95.9
2 2 56.8 4.9 6.5 4.2 3.8 2 12.60 8 45 18 8.0 5 49 16 19 3 2 9.6 93.5
3 3 54.8 0.2 5.4 6.1 6.7 2 8.29 12 41 21 14.7 2 39 17 15 1 2 1.7 95.5

✏️ Consigna 1

Uní las dos tablas para tener, en una sola fila por empleado, su salario y su bienestar.

a) salario tiene tres años por persona. Quedate con 2025.

b) Unilas por la columna que comparten. Las dos tienen empleado_id.

salario_2025 <- salario |>
  # TODO: completá el año que nos interesa
  filter(anio == ___) |>
  select(empleado_id, salario_mensual)

datos <- clima |>
  # TODO: completá la columna por la que se unen las dos tablas
  inner_join(salario_2025, by = "___") |>
  mutate(salario = salario_mensual / 1e6)   # en millones, para leerlo cómodo

cat(dim(datos), "\n")
head(datos[, c("empleado_id", "salario", "bienestar_laboral")], 3)
correlacion <- cor(datos$salario, datos$bienestar_laboral)
cat(sprintf("correlación salario / bienestar: %.3f\n", correlacion))
summary(datos[, c("salario", "bienestar_laboral")])
correlación salario / bienestar: 0.572
    salario      bienestar_laboral
 Min.   :0.935   Min.   :11.40    
 1st Qu.:1.191   1st Qu.:49.90    
 Median :1.264   Median :57.20    
 Mean   :1.258   Mean   :57.36    
 3rd Qu.:1.323   3rd Qu.:64.50    
 Max.   :1.522   Max.   :97.10

El punto de partida: el scatter

# Las dos variables sueltas: se usan en casi todas las celdas de acá en adelante.
x <- datos$salario
y <- datos$bienestar_laboral

nube <- ggplot(datos, aes(salario, bienestar_laboral)) +
  geom_point(alpha = 0.35, colour = "#33404a", size = 1.4) +
  labs(x = "Salario mensual (millones de $)", y = "Bienestar laboral (0-100)") +
  ylim(0, 100) +
  theme_minimal(base_size = 12)

nube + ggtitle("Nimbus: 600 empleados")

2. Una recta a ojo

Una recta son dos números: la ordenada al origen β₀ y la pendiente β₁.

\[\hat{y} = \beta_0 + \beta_1 x\]

donde - \(\hat{y}\) es el bienestar que la recta predice, - \(x\) es el salario, - \(\beta_0\) es dónde cruza el eje vertical, - \(\beta_1\) es cuánto sube el bienestar por cada millón más de salario.

En R no hay sliders como en Colab con Python, así que probamos tres rectas a la vez y comparamos. El RSS es la suma de los residuos al cuadrado, y es lo que dice cuál es mejor:

\[\text{RSS} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]

✏️ Consigna 2

Programá el RSS: calculá el residuo de cada persona y sumá todos esos residuos al cuadrado.

rss <- function(b0, b1) {
  # TODO: el residuo es el valor real menos el que predice la recta b0 + b1*x
  residuos <- y - (___ + ___ * x)
  # TODO: ¿a qué potencia hay que elevarlos antes de sumar?
  sum(residuos^___)
}

rss_mala <- rss(25, 26)
cat(sprintf("RSS de la recta (25,0 · 26,0): %s\n", format(round(rss_mala), big.mark = ".")))
# Tres candidatas, para ver que el RSS ordena lo que el ojo intuye.
candidatas <- data.frame(
  nombre = c("muy plana", "razonable", "muy empinada"),
  b0 = c(45, 0, -70),
  b1 = c(10, 45, 100)
)
candidatas$rss <- mapply(rss, candidatas$b0, candidatas$b1)
candidatas$rss <- round(candidatas$rss)
candidatas
A data.frame: 3 × 4
nombre b0 b1 rss
<chr> <dbl> <dbl> <dbl>
muy plana 45 10 69119
razonable 0 45 54438
muy empinada -70 100 56839 
nube +
  geom_abline(data = candidatas, aes(intercept = b0, slope = b1, colour = nombre), linewidth = 1.1) +
  scale_colour_manual(values = c("muy plana" = "#C8622A", "razonable" = "#00529B",
                                 "muy empinada" = "#1F7A4D"), name = NULL) +
  ggtitle("Tres rectas elegidas a ojo") +
  theme(legend.position = "bottom")

3. Mínimos cuadrados, sin lm()

Para este problema hay una fórmula exacta:

\[\hat{\beta}_1 = \frac{\sum_i (x_i - \bar{x})(y_i - \bar{y})}{\sum_i (x_i - \bar{x})^2} \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\]

donde \(\bar{x}\) y \(\bar{y}\) son los promedios. La segunda dice algo lindo: la recta pasa siempre por el punto de los dos promedios.

✏️ Consigna 3

Programá las dos fórmulas de arriba. Las dos usan los promedios mx y my, que ya están calculados en la primera línea.

Más adelante comparamos este resultado contra lm(): si ahí coinciden, programaste bien la fórmula.

minimos_cuadrados <- function(x, y) {
  mx <- mean(x); my <- mean(y)
  # TODO: arriba, cuánto se aparta cada uno de SU promedio; abajo, sólo el de x
  b1 <- sum((x - mx) * (y - ___)) / sum((x - ___)^2)
  # TODO: la recta pasa por (mx, my). Despejá b0 de my = b0 + b1 * mx
  b0 <- my - ___ * mx
  c(b0 = b0, b1 = b1)
}

ols <- minimos_cuadrados(x, y)
b0_ols <- unname(ols["b0"]); b1_ols <- unname(ols["b1"])
cat(sprintf("β₀ = %.2f\nβ₁ = %.2f\n", b0_ols, b1_ols))
# Si nos movemos un poco en cualquier dirección, el RSS empeora.
set.seed(42)
rss_ols <- rss(b0_ols, b1_ols)

for (db in list(c(-5, 0), c(5, 0), c(0, -5), c(0, 5))) {
  cat(sprintf("  β₀%+.0f  β₁%+.0f   RSS = %s\n", db[1], db[2],
              format(round(rss(b0_ols + db[1], b1_ols + db[2])), big.mark = ".")))
}
cat(sprintf("  la de mínimos cuadrados   RSS = %s   <- la más chica\n",
            format(round(rss_ols), big.mark = ".")))

# La recta pasa por el punto de los dos promedios.
cat(sprintf("\npredicción en el salario promedio: %.2f\n", b0_ols + b1_ols * mean(x)))
cat(sprintf("promedio real de bienestar:        %.2f\n", mean(y)))
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
  β₀-5  β₁+0   RSS = 65.874
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
  β₀+5  β₁+0   RSS = 65.874
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
  β₀+0  β₁-5   RSS = 74.722
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
  β₀+0  β₁+5   RSS = 74.722
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
  la de mínimos cuadrados   RSS = 50.874   <- la más chica

predicción en el salario promedio: 57.36
promedio real de bienestar:        57.36
nube +
  geom_abline(intercept = 0, slope = 45, colour = "#C8622A", linewidth = 1.1, linetype = "dashed") +
  geom_abline(intercept = b0_ols, slope = b1_ols, colour = "#00529B", linewidth = 1.3) +
  ggtitle(sprintf("A ojo (RSS = %s) contra mínimos cuadrados (RSS = %s)",
                  format(round(rss(0, 45)), big.mark = "."),
                  format(round(rss_ols), big.mark = ".")))
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”
Warning message in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
“'big.mark' and 'decimal.mark' are both '.', which could be confusing”

4. Medir el ajuste: RSE y R²

\[\text{RSE} = \sqrt{\frac{\text{RSS}}{n - 2}} \qquad R^2 = 1 - \frac{\text{RSS}}{\text{TSS}} \qquad \text{TSS} = \sum_i (y_i - \bar{y})^2\]

donde el \(-2\) del RSE es porque estimamos dos parámetros, y \(R^2 = 0\) quiere decir que la recta no le ganó a decir el promedio.

✏️ Consigna 4

Calculá las dos medidas. El TSS ya está hecho. Ojo con el denominador del RSE: no es n, es n menos la cantidad de parámetros que estimamos.

n <- length(y)
tss <- sum((y - mean(y))^2)

# TODO: ¿cuántos parámetros estimamos en una regresión simple?
rse <- sqrt(rss_ols / (n - ___))
# TODO: el R² compara el error de nuestra recta contra el del modelo del promedio
r2 <- 1 - ___ / tss

cat(sprintf("TSS = %s\nRSE = %.2f puntos de bienestar\nR²  = %.3f\n",
            format(round(tss), big.mark = "."), rse, r2))

Lo mismo con lm()

En R el ajuste es una línea. Lo importante es verificar que da exactamente lo mismo.

Un detalle de vocabulario: lo que summary() llama Residual standard error es nuestro RSE, y Multiple R-squared es el R².

modelo <- lm(bienestar_laboral ~ salario, data = datos)
resumen <- summary(modelo)

cat(sprintf("lm()   β₀ = %.4f   β₁ = %.4f   R² = %.4f   RSE = %.4f\n",
            coef(modelo)[[1]], coef(modelo)[[2]], resumen$r.squared, resumen$sigma))
cat(sprintf("a mano β₀ = %.4f   β₁ = %.4f   R² = %.4f   RSE = %.4f\n",
            b0_ols, b1_ols, r2, rse))
resumen
lm()   β₀ = -31.0448   β₁ = 70.2997   R² = 0.3274   RSE = 9.2235
a mano β₀ = -31.0448   β₁ = 70.2997   R² = 0.3274   RSE = 9.2235

Call:
lm(formula = bienestar_laboral ~ salario, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.187  -5.602   0.460   5.385  32.050 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -31.045      5.196  -5.975 3.95e-09 ***
salario       70.300      4.121  17.061  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.224 on 598 degrees of freedom
Multiple R-squared:  0.3274,    Adjusted R-squared:  0.3263 
F-statistic: 291.1 on 1 and 598 DF,  p-value: < 2.2e-16

En la salida de summary() fijate en la columna Pr(>|t|) del salario: es el p-valor del que hablamos en la teoría. Un número minúsculo quiere decir que el cero no es una explicación razonable de la pendiente que estamos viendo.

5. Varias variables

La encuesta tiene 19 predictores además del salario. Estos son los candidatos, cada uno solo, con su R²:

CANDIDATAS <- c("salario", "apoyo_equipo", "reconocimiento", "autonomia",
                "horas_extra_semana", "dias_home_office", "bono_anual_pct",
                "reuniones_semana", "mensajes_chat_dia", "puntualidad_pct",
                "horas_capacitacion", "distancia_oficina_km")

r2_de <- function(columnas, tabla = datos) {
  # R² del modelo que usa esas columnas. Es la función que vas a usar todo el rato.
  summary(lm(reformulate(columnas, response = "bienestar_laboral"), data = tabla))$r.squared
}

solas <- c()
for (columna in CANDIDATAS) {
  solas[columna] <- r2_de(columna)
}
round(sort(solas, decreasing = TRUE), 3)
salario
0.327
bono_anual_pct
0.306
apoyo_equipo
0.16
horas_extra_semana
0.082
reconocimiento
0.076
autonomia
0.007
horas_capacitacion
0.006
reuniones_semana
0.004
dias_home_office
0.002
puntualidad_pct
0.001
mensajes_chat_dia
0
distancia_oficina_km
0

Hay tres que explican algo y el resto está cerca de cero. La pregunta obvia: si junto las tres buenas, ¿el R² es la suma de los tres?

✏️ Consigna 5

Compará las dos cosas: sumar los tres R² por separado contra ajustar un solo modelo con las tres variables juntas.

Anotá tu predicción antes de correr la celda: ¿va a dar más, menos o igual?

tres <- c("salario", "apoyo_equipo", "reconocimiento")

suma_individuales <- sum(solas[tres])
# TODO: ajustá UN modelo con las tres variables a la vez
juntas <- r2_de(___)

cat(sprintf("sumando los tres R2 por separado:   %.3f\n", suma_individuales))
cat(sprintf("el modelo con las tres juntas:      %.3f\n", juntas))
cat(sprintf("diferencia (información repetida):  %.3f\n", suma_individuales - juntas))

No se suman: la información se superpone.

Un coeficiente no se lee solo

solo_bono <- lm(bienestar_laboral ~ bono_anual_pct, data = datos)
con_salario <- lm(bienestar_laboral ~ bono_anual_pct + salario, data = datos)
r_bono_salario <- cor(datos$bono_anual_pct, datos$salario)

cat(sprintf("beta del bono, solo:           %+.3f\n", coef(solo_bono)[["bono_anual_pct"]]))
cat(sprintf("beta del bono, con el salario: %+.3f\n", coef(con_salario)[["bono_anual_pct"]]))
cat(sprintf("correlación bono / salario:     %.3f\n", r_bono_salario))
cat(sprintf("R2 solo bono: %.3f   con las dos: %.3f\n",
            summary(solo_bono)$r.squared, summary(con_salario)$r.squared))
beta del bono, solo:           +1.845
beta del bono, con el salario: +0.010
correlación bono / salario:     0.966
R2 solo bono: 0.306   con las dos: 0.327

El bono parecía un predictor fuerte y resultó ser el salario escrito otra vez: en Nimbus el bono es un porcentaje del sueldo.

Ahora armá tu modelo

Esta es la parte en la que trabajás vos. Tenés CANDIDATAS y r2_de().

# Un ayudante para ir tanteando rápido.
probar <- function(columnas) {
  valor <- r2_de(columnas)
  cat(sprintf("R2 = %.3f   con %d variable(s): %s\n",
              valor, length(columnas), paste(columnas, collapse = ", ")))
  invisible(valor)
}

probar(c("salario"))
probar(c("salario", "apoyo_equipo"))
probar(c("salario", "apoyo_equipo", "reconocimiento"))
probar(c("salario", "apoyo_equipo", "reconocimiento", "puntualidad_pct"))
cat("\nVariables disponibles:\n")
cat(paste(CANDIDATAS, collapse = ", "), "\n")
R2 = 0.327   con 1 variable(s): salario
R2 = 0.501   con 2 variable(s): salario, apoyo_equipo
R2 = 0.522   con 3 variable(s): salario, apoyo_equipo, reconocimiento
R2 = 0.523   con 4 variable(s): salario, apoyo_equipo, reconocimiento, puntualidad_pct

Variables disponibles:
salario, apoyo_equipo, reconocimiento, autonomia, horas_extra_semana, dias_home_office, bono_anual_pct, reuniones_semana, mensajes_chat_dia, puntualidad_pct, horas_capacitacion, distancia_oficina_km

Mirá las cuatro pruebas antes de seguir. Agregar apoyo_equipo y reconocimiento sube el R². Agregar puntualidad_pct, que es ruido puro, lo sube también, aunque poquísimo. Guardate eso: en un rato va a ser el centro de la clase.

✏️ Consigna 6

Armá tu propio modelo con al menos tres variables de CANDIDATAS, usando probar() para ir tanteando. Cuando tengas una combinación que te convenza, ponela en mis_variables.

No busques la mejor de todas: buscá una que te parezca razonable.

# TODO: poné acá las variables que elegiste. Al menos tres, todas de CANDIDATAS.
mis_variables <- c("___", "___", "___")

r2_mio <- r2_de(mis_variables)
probar(mis_variables)

6. El número honesto: entrenamiento y testeo

Todo lo que hiciste hasta acá lo mediste sobre las mismas 600 personas con las que armaste el modelo. Es como corregir un examen con el machete a la vista.

✏️ Consigna 7

Partí la tabla dejando el 30% para testeo, y evaluá tu modelo en los dos lados.

set.seed(42) es para que a todos les dé lo mismo y podamos comparar en clase.

set.seed(42)
# TODO: qué proporción va a ENTRENAMIENTO (el resto queda para testeo)
indices <- sample(nrow(datos), size = round(___ * nrow(datos)))
entrenamiento <- datos[indices, ]
testeo <- datos[-indices, ]

r2_fuera <- function(modelo, evaluacion) {
  pred <- predict(modelo, newdata = evaluacion)
  1 - sum((evaluacion$bienestar_laboral - pred)^2) /
    sum((evaluacion$bienestar_laboral - mean(entrenamiento$bienestar_laboral))^2)
}

mio <- lm(reformulate(mis_variables, response = "bienestar_laboral"), data = entrenamiento)
r2_train_mio <- summary(mio)$r.squared
# TODO: el número honesto se mide en el conjunto que el modelo NO vio
r2_test_mio <- r2_fuera(mio, ___)

cat(sprintf("tu modelo: %s\n", paste(mis_variables, collapse = ", ")))
cat(sprintf("  filas de entrenamiento: %d   filas de testeo: %d\n", nrow(entrenamiento), nrow(testeo)))
cat(sprintf("  R2 con TODA la data (lo de recién): %.3f\n", r2_mio))
cat(sprintf("  R2 en entrenamiento:                %.3f\n", r2_train_mio))
cat(sprintf("  R2 en testeo:                       %.3f\n", r2_test_mio))
cat(sprintf("  hueco entrenamiento - testeo:       %+.3f\n", r2_train_mio - r2_test_mio))

¿Se podía elegir mejor?

Comparemos tres modelos sobre la misma partición: sólo el salario, el tuyo, y uno con todas las variables.

Mirá las dos columnas por separado. Si agregar variables siempre mejorara, el de todas tendría que ganar en las dos. Prestá atención a la última columna, el hueco entre entrenamiento y testeo: es la medida de cuánto se está engañando cada modelo a sí mismo.

TODAS <- setdiff(names(datos), c("empleado_id", "bienestar_laboral", "salario_mensual"))

evaluar <- function(columnas) {
  m <- lm(reformulate(columnas, response = "bienestar_laboral"), data = entrenamiento)
  c(train = summary(m)$r.squared, test = r2_fuera(m, testeo))
}

nombres <- c("sólo el salario", "el tuyo", "TODAS")
listas_columnas <- list(c("salario"), mis_variables, TODAS)

comparacion <- data.frame()
for (i in seq_along(nombres)) {
  r <- evaluar(listas_columnas[[i]])
  fila <- data.frame(modelo = nombres[i], variables = length(listas_columnas[[i]]),
                     train = r[["train"]], test = r[["test"]], hueco = r[["train"]] - r[["test"]])
  comparacion <- rbind(comparacion, fila)
}

comparacion$train <- round(comparacion$train, 3)
comparacion$test <- round(comparacion$test, 3)
comparacion$hueco <- round(comparacion$hueco, 3)
comparacion
A data.frame: 3 × 5
modelo variables train test hueco
<chr> <int> <dbl> <dbl> <dbl>
sólo el salario 1 0.360 0.241 0.119
el tuyo 4 0.542 0.526 0.016
TODAS 19 0.630 0.589 0.042

Ahí está la clase entera en una tabla. Con todas las variables el R² de entrenamiento es el más alto de los tres, y el de testeo no. El hueco entre las dos columnas es el sobreajuste.

La curva completa

Para ver dónde estaba el óptimo, armamos una escalera: empezamos con la variable que más explicaba sola, y le vamos sumando las demás en el orden en que aparecían en esa lista, de la que más explicaba a la que menos. En cada escalón medimos de los dos lados.

Es a mano y a propósito. Hay 4.095 combinaciones posibles de estas doce variables, y no hace falta recorrerlas para ver lo que queremos ver.

# `solas` tiene el R2 de cada variable sola. La ordenamos de mayor a menor y las
# vamos agregando en ese orden, una por escalón.
orden <- names(sort(solas, decreasing = TRUE))

historia <- data.frame()
for (k in 1:length(orden)) {
  r <- evaluar(orden[1:k])                       # las primeras k de la lista
  fila <- data.frame(n_variables = k, agregada = orden[k],
                     r2_entrenamiento = r[["train"]], r2_testeo = r[["test"]])
  historia <- rbind(historia, fila)
}

mejor <- historia[which.max(historia$r2_testeo), ]

cat(sprintf("mejor en testeo: %d variables, R2 = %.3f\n", mejor$n_variables, mejor$r2_testeo))
cat(sprintf("tu modelo:       %d variables, R2 = %.3f\n", length(mis_variables), r2_test_mio))

# round() no acepta un data.frame con columnas de texto, así que redondeamos
# solo las dos numéricas.
historia$r2_entrenamiento <- round(historia$r2_entrenamiento, 3)
historia$r2_testeo <- round(historia$r2_testeo, 3)
historia
mejor en testeo: 6 variables, R2 = 0.604
tu modelo:       4 variables, R2 = 0.526
A data.frame: 12 × 4
n_variables agregada r2_entrenamiento r2_testeo
<int> <chr> <dbl> <dbl>
1 salario 0.360 0.241
2 bono_anual_pct 0.360 0.239
3 apoyo_equipo 0.522 0.447
4 horas_extra_semana 0.605 0.523
5 reconocimiento 0.614 0.566
6 autonomia 0.621 0.604
7 horas_capacitacion 0.621 0.603
8 reuniones_semana 0.622 0.599
9 dias_home_office 0.625 0.596
10 puntualidad_pct 0.626 0.593
11 mensajes_chat_dia 0.627 0.590
12 distancia_oficina_km 0.627 0.591 
largo <- rbind(
  data.frame(n = historia$n_variables, r2 = historia$r2_entrenamiento, conjunto = "entrenamiento"),
  data.frame(n = historia$n_variables, r2 = historia$r2_testeo, conjunto = "testeo")
)

ggplot(largo, aes(n, r2, colour = conjunto)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.8) +
  geom_vline(xintercept = mejor$n_variables, linetype = "dashed", colour = "#33404a") +
  annotate("point", x = length(mis_variables), y = r2_test_mio,
           shape = 8, size = 4, colour = "#1F7A4D", stroke = 1.2) +
  scale_colour_manual(values = c(entrenamiento = "#1F6FB4", testeo = "#C8622A")) +
  labs(x = "Variables en el modelo", y = "R2", colour = NULL,
       title = "El entrenamiento nunca baja. El testeo sí.",
       subtitle = "La estrella verde es tu modelo, medido en testeo") +
  theme_minimal(base_size = 12)

La curva azul sube siempre; la naranja sube, hace un máximo y después se cae. Ese máximo es la cantidad de variables que conviene usar, y no es la mayor. La estrella verde es dónde quedó tu modelo.

Fijate en el escalón del bono: agregarlo no mueve ninguno de los dos R². Es la misma historia de hace un rato, ahora medida. Una variable puede tener un R² alto ella sola y aportar exactamente cero al lado de las demás.

Lo que queda: penalizar en vez de descartar

glmnet ajusta ridge (alpha = 0) y lasso (alpha = 1). Estandariza por defecto, que es lo que hay que hacer siempre antes de penalizar.

Los números no coinciden dígito a dígito con los de Python, y los lambda tampoco. glmnet parametriza la penalización de otra manera que scikit-learn: acá lambda = 10 para ridge y lambda = 1 para lasso hacen el mismo trabajo que alpha = 10 y alpha = 0,5 allá. Lo que sí coincide es el orden: con pocos datos mínimos cuadrados se desarma y las dos penalizadas aguantan.

✏️ Consigna 8

Comprobá con los datos lo que dijimos en la teoría. Entrená los tres modelos con sólo 30 personas y las 19 variables, y comparalos sobre el resto.

En glmnet, alpha = 0 es ridge y alpha = 1 es lasso. No confundas ese alpha con el alpha de Python, que ahí es el lambda.

library(glmnet)

set.seed(11)
# TODO: cuántas personas usamos para entrenar
chicos <- sample(nrow(datos), size = ___)
chico <- datos[chicos, ]
resto <- datos[-chicos, ]

X_chico <- as.matrix(chico[, TODAS]); X_resto <- as.matrix(resto[, TODAS])
y_chico <- chico$bienestar_laboral;   y_resto <- resto$bienestar_laboral

r2_resto <- function(pred) 1 - sum((y_resto - pred)^2) / sum((y_resto - mean(y_chico))^2)

ols_chico <- lm(reformulate(TODAS, response = "bienestar_laboral"), data = chico)
# TODO: alpha = 0 es ridge, alpha = 1 es lasso
ridge <- glmnet(X_chico, y_chico, alpha = ___, lambda = 10)
lasso <- glmnet(X_chico, y_chico, alpha = ___, lambda = 1)

resultados <- c(
  "mínimos cuadrados" = r2_resto(predict(ols_chico, newdata = resto)),
  "Ridge"             = r2_resto(as.numeric(predict(ridge, newx = X_resto))),
  "Lasso"             = r2_resto(as.numeric(predict(lasso, newx = X_resto)))
)

for (nombre in names(resultados)) {
  cat(sprintf("%18s: R2 de testeo = %+.3f\n", nombre, resultados[[nombre]]))
}

Con 30 personas y 19 variables, mínimos cuadrados da un R² negativo: predice peor que decir el promedio y no pensar.

Hoja de referencia

Qué Fórmula En R
Residuo \(e_i = y_i - \hat{y}_i\) residuals(modelo)
RSS \(\sum_i e_i^2\) sum(residuals(modelo)^2)
RSE \(\sqrt{\text{RSS}/(n-2)}\) summary(modelo)$sigma
R² \(1 - \text{RSS}/\text{TSS}\) summary(modelo)$r.squared
Ajustar lm(y ~ x1 + x2, data = datos)
Predecir predict(modelo, newdata = testeo)
Ridge RSS \(+ \lambda \sum \beta_j^2\) glmnet(X, y, alpha = 0)
Lasso RSS \(+ \lambda \sum \lvert\beta_j\rvert\) glmnet(X, y, alpha = 1)

Cuidado con var() y sd() en R: dividen por \(n-1\). En numpy el default es \(n\). Acá no molesta porque todo se calculó con sumas explícitas, pero si comparás resultados entre los dos lenguajes es la primera cosa a revisar.

Para el trabajo práctico

Sobre el dataset de tu grupo:

  1. Elegí una variable respuesta numérica.
  2. Ajustá la regresión simple con el predictor más prometedor y reportá el R² de testeo.
  3. Agregá predictores de a uno y armá la curva de R² de entrenamiento contra R² de testeo, como la figura de la sección 6. Marcá dónde se separan.
  4. Una línea de conclusión: ¿el mejor modelo es el que más variables tiene? ¿Por qué?