2.6 - Regresión lineal múltiple: el modelo de mascarillas de cinco variables¶
En el capítulo 2, Marta extiende la recta de la regresión lineal a un plano y, con más variables, a un hiperplano: $\hat{y} = w_1 x_1 + w_2 x_2 + \dots + w_n x_n + b$. Carlos aporta dieciocho meses de registros del hospital, Carmen propone las variables clínicas que importan y el equipo ajusta un modelo que predice las mascarillas necesarias cada día a partir de cinco variables: pacientes, temperatura, urgencias, alertas sanitarias y fines de semana.
Los otros dos cuadernos del capítulo usan una sola variable. Aquí implementamos el modelo con varias, como lo describe el libro. Para poder comprobar el resultado, generamos los datos a partir de la ecuación que obtuvo el equipo en el libro, más algo de ruido, y vemos si cada método la recupera:
- Importar las bibliotecas.
- De la recta al hiperplano: el ejemplo de Marta.
- Dieciocho meses de datos.
- Meses de entrenamiento y de prueba.
- El descenso del gradiente con las ecuaciones de actualización de Luis.
- Escalar las características.
- La solución exacta: mínimos cuadrados y scikit-learn.
- ¿Cómo de bueno es el modelo?
- Cuando los datos tienen errores: ¿ECM o EAM?
- La previsión para la semana que viene.
Paso 1: Importar las bibliotecas¶
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, QuantileRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
np.set_printoptions(precision=3, suppress=True)
Paso 2: De la recta al hiperplano¶
Con varias variables, la predicción es una suma ponderada: cada variable $x_i$ se multiplica por su peso $w_i$ y se suma el sesgo $b$. Escrito con vectores, $\hat{y} = \mathbf{w} \cdot \mathbf{x} + b$, que es exactamente lo que calcula el operador @ de NumPy. El ejemplo de Marta en el libro usa tres variables, pacientes, temperatura y nivel de alerta sanitaria, con $\hat{y} = 0.6\,x_1 - 0.4\,x_2 + 20\,x_3 + 50$, para un día con 300 pacientes, 8 °C y sin alerta:
w = np.array([0.6, -0.4, 20.0]) # mascarillas por paciente, por grado, por nivel de alerta
b = 50.0 # la base diaria
x = np.array([300, 8, 0]) # 300 pacientes, 8 grados, sin alerta
y_hat = w @ x + b
print(f"Mascarillas previstas: {y_hat:.1f}")
print("Término a término:", w * x, "+", b)
Mascarillas previstas: 226.8 Término a término: [180. -3.2 0. ] + 50.0
226,8 mascarillas: 180 por los 300 pacientes, −3,2 por la temperatura, nada por la alerta, más la base de 50. Como dice Marta, el decimal no significa una fracción de mascarilla; es una estimación, y en la práctica el centro prepararía unas 227. Una predicción con cinco variables es la misma operación con vectores más largos.
Paso 3: Dieciocho meses de datos¶
Los registros de Carlos cubren dieciocho meses. Los simulamos, con una fila por día, con las variables del modelo que ajusta el equipo en el libro: los pacientes atendidos (menos en fin de semana), la temperatura (más fría en invierno), las urgencias (más en invierno), tres periodos de alerta sanitaria de dos semanas y si el día es de fin de semana. Las mascarillas usadas se generan con la ecuación del equipo,
$$\text{mascarillas} = 42.3 + 0.05\,\text{pacientes} - 0.37\,\text{temperatura} + 3.2\,\text{urgencias} + 18.7\,\text{alerta} + 4.1\,\text{fin de semana}$$
más un ruido aleatorio con una desviación típica de 4 mascarillas. Los modelos solo verán los datos, nunca esta ecuación, así que después podremos comprobar cuánto se acerca cada uno a ella.
rng = np.random.default_rng(7)
dates = pd.date_range('2024-01-01', '2025-06-30', freq='D')
n = len(dates)
day_of_year = dates.dayofyear.values
winter = np.cos(2 * np.pi * (day_of_year - 20) / 365) # 1 a mediados de enero, -1 a mediados de julio
fin_de_semana = (dates.dayofweek >= 5).astype(int)
temperatura = 16 - 8 * winter + rng.normal(0, 2.5, n)
pacientes = np.where(fin_de_semana == 1, rng.normal(245, 18, n), rng.normal(295, 18, n)).round()
urgencias = rng.poisson(np.where(winter > 0.5, 13, 10))
alerta = np.zeros(n, dtype=int)
for start in ['2024-01-15', '2024-12-20', '2025-02-03']:
first = dates.get_loc(pd.Timestamp(start))
alerta[first:first + 14] = 1
book_w = np.array([0.05, -0.37, 3.2, 18.7, 4.1])
book_b = 42.3
X = np.c_[pacientes, temperatura, urgencias, alerta, fin_de_semana]
mascarillas = X @ book_w + book_b + rng.normal(0, 4, n)
features = ['pacientes', 'temperatura', 'urgencias', 'alerta', 'fin_de_semana']
data = pd.DataFrame(X, index=dates, columns=features).assign(mascarillas=mascarillas)
print(data.round(1).head())
print(f"\n{n} días, {alerta.sum()} de ellos con alerta sanitaria")
pacientes temperatura urgencias alerta fin_de_semana \
2024-01-01 309.0 8.4 9.0 0.0 0.0
2024-01-02 301.0 9.1 12.0 0.0 0.0
2024-01-03 283.0 7.7 15.0 0.0 0.0
2024-01-04 332.0 6.1 16.0 0.0 0.0
2024-01-05 337.0 7.1 13.0 0.0 0.0
mascarillas
2024-01-01 91.6
2024-01-02 83.4
2024-01-03 101.9
2024-01-04 109.6
2024-01-05 93.9
547 días, 42 de ellos con alerta sanitaria
fig, axes = plt.subplots(1, 3, figsize=(16, 4))
axes[0].plot(data.index, data['mascarillas'], lw=0.8)
axes[0].set_title('Mascarillas usadas cada día')
axes[0].tick_params(axis='x', rotation=30)
axes[1].scatter(data['pacientes'], data['mascarillas'], s=8, c=data['fin_de_semana'], cmap='coolwarm')
axes[1].set_xlabel('pacientes'); axes[1].set_title('Mascarillas y pacientes (rojo: fin de semana)')
axes[2].scatter(data['temperatura'], data['mascarillas'], s=8, c=data['alerta'], cmap='coolwarm')
axes[2].set_xlabel('temperatura (°C)'); axes[2].set_title('Mascarillas y temperatura (rojo: alerta)')
plt.tight_layout()
plt.show()
547 días, 42 de ellos con alerta sanitaria. Las mascarillas suben en los dos inviernos, cuando hace más frío y hay más urgencias, y se disparan durante las alertas (en rojo en la gráfica de la derecha). El número de pacientes, por sí solo, dice poco: su peso en la ecuación del libro es pequeño, y los fines de semana, con menos pacientes, siguen usando las mismas mascarillas. Por eso el modelo necesita todas las variables a la vez.
Paso 4: Meses de entrenamiento y de prueba¶
La regla de Alma en el capítulo 2 es no evaluar nunca un modelo con los datos con los que se entrenó. Con datos ordenados en el tiempo, la separación honesta es por fechas: entrenamos con los quince primeros meses y guardamos los tres últimos, que el modelo no verá nunca, para probarlo. Así es también como lo usará el hospital: entrenado con el pasado, preguntado por el futuro.
train = dates < pd.Timestamp('2025-04-01')
X_train, y_train = X[train], mascarillas[train]
X_test, y_test = X[~train], mascarillas[~train]
m = len(y_train)
print(f"Entrenamiento: {m} días (del {dates[train][0].date()} al {dates[train][-1].date()})")
print(f"Prueba: {len(y_test)} días (del {dates[~train][0].date()} al {dates[~train][-1].date()})")
Entrenamiento: 456 días (del 2024-01-01 al 2025-03-31) Prueba: 91 días (del 2025-04-01 al 2025-06-30)
Paso 5: El descenso del gradiente con las ecuaciones de actualización de Luis¶
Las notas de Luis dan la actualización de cada peso y del sesgo para la regresión lineal con el coste ECM:
$$w_j := w_j - \alpha \, \frac{2}{m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})\, x_j^{(i)} \qquad b := b - \alpha \, \frac{2}{m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})$$
Con una matriz $X$ de $m$ filas y una columna por variable, las cinco sumas de los pesos son un único producto, $X^\top (\hat{y} - y)$. Primero las aplicamos a las variables tal como vienen. El descenso del gradiente diverge si la tasa de aprendizaje es demasiado grande, y con estas variables se puede calcular la mayor tasa que todavía es estable: depende de la variable más grande, el número de pacientes, que está en los cientos.
def gradient_descent(X, y, learning_rate, iterations):
m = len(y)
w, b = np.zeros(X.shape[1]), 0.0
for _ in range(iterations):
error = X @ w + b - y # y_hat - y para cada día
w = w - learning_rate * (2 / m) * (X.T @ error) # la actualización de Luis, todos los pesos a la vez
b = b - learning_rate * (2 / m) * error.sum()
return w, b
def mse(X, y, w, b):
return np.mean((X @ w + b - y) ** 2)
A = np.c_[X_train, np.ones(m)]
largest = np.linalg.eigvalsh((2 / m) * A.T @ A).max()
print(f"Mayor tasa de aprendizaje estable con las variables sin escalar: {2 / largest:.2e}")
w_raw, b_raw = gradient_descent(X_train, y_train, learning_rate=0.95 * 2 / largest, iterations=5000)
print("\nTras 5000 iteraciones al 95% de esa tasa:")
print(" pesos:", w_raw, " sesgo:", round(b_raw, 2))
print(f" ECM de entrenamiento: {mse(X_train, y_train, w_raw, b_raw):.2f}")
Mayor tasa de aprendizaje estable con las variables sin escalar: 1.24e-05 Tras 5000 iteraciones al 95% de esa tasa: pesos: [ 0.206 -0.432 3.34 0.201 0.376] sesgo: 0.09 ECM de entrenamiento: 76.17
La mayor tasa de aprendizaje estable es de solo 1,24e-05, porque los pacientes están en los cientos mientras que la alerta y el fin de semana valen 0 o 1. A esa tasa, tras 5000 iteraciones, el modelo sigue lejos de los datos: el peso de la alerta es 0,201 en lugar de unos 18, el sesgo es 0,09 en lugar de unos 44 y el ECM de entrenamiento es 76,17. Cada paso mueve mucho el peso de los pacientes y casi nada el resto, así que el descenso se arrastra por un valle largo y estrecho. Es el problema que resuelve el cuaderno desde cero de este capítulo escalando su única variable, ahora con cinco.
Paso 6: Escalar las características¶
Las notas de Luis advierten de que el descenso del gradiente es sensible a la escala de las variables, y recomiendan normalizar los datos. Estandarizamos cada variable: le restamos su media y la dividimos por su desviación típica, las dos calculadas solo con los días de entrenamiento. Así todas las variables tienen media 0 y desviación típica 1, la misma tasa de aprendizaje sirve para todos los pesos y el descenso puede dar pasos grandes. Al final devolvemos los pesos y el sesgo a las unidades originales, para que cada peso vuelva a ser "mascarillas por paciente", "mascarillas por grado", etc.:
$$w_j = \frac{w'_j}{\sigma_j} \qquad b = b' - \sum_j \frac{w'_j\,\mu_j}{\sigma_j}$$
mu, sigma = X_train.mean(axis=0), X_train.std(axis=0)
Z_train = (X_train - mu) / sigma
for iterations in [10, 50, 200, 1000]:
w_s, b_s = gradient_descent(Z_train, y_train, learning_rate=0.1, iterations=iterations)
print(f"{iterations:>5} iteraciones: ECM de entrenamiento {mse(Z_train, y_train, w_s, b_s):8.3f}")
w_gd = w_s / sigma # de vuelta a mascarillas por paciente, por grado, ...
b_gd = b_s - np.sum(w_s * mu / sigma)
print("\nPesos en unidades originales:", w_gd, " sesgo:", round(b_gd, 2))
10 iteraciones: ECM de entrenamiento 112.756 50 iteraciones: ECM de entrenamiento 15.246 200 iteraciones: ECM de entrenamiento 15.238 1000 iteraciones: ECM de entrenamiento 15.238 Pesos en unidades originales: [ 0.044 -0.374 3.207 18.808 3.363] sesgo: 44.3
Con las variables estandarizadas, la misma función llega a un ECM de entrenamiento de 15,246 en 50 iteraciones y de 15,238 a partir de 200, cerca de la varianza del ruido que añadimos (16). Devuelto a las unidades originales, el modelo necesita 0,044 mascarillas por paciente, −0,374 por grado, 3,207 por urgencia, 18,808 más en un día de alerta y 3,363 más en fin de semana, sobre una base de 44,3.
Paso 7: La solución exacta: mínimos cuadrados y scikit-learn¶
Para la regresión lineal con el coste ECM existe también una solución directa, el método de los mínimos cuadrados que menciona Marta: los pesos que minimizan el error al cuadrado resuelven un sistema de ecuaciones lineales, las ecuaciones normales $(A^\top A)\,\theta = A^\top y$, donde $A$ es $X$ con una columna extra de unos para el sesgo. NumPy lo resuelve con una sola llamada, y el LinearRegression de scikit-learn, que usa el equipo en el libro, hace lo mismo dentro de fit:
theta = np.linalg.lstsq(A, y_train, rcond=None)[0]
model = LinearRegression().fit(X_train, y_train)
table = pd.DataFrame({'ecuación del libro': np.r_[book_w, book_b],
'descenso del gradiente': np.r_[w_gd, b_gd],
'mínimos cuadrados': theta,
'scikit-learn': np.r_[model.coef_, model.intercept_]},
index=features + ['sesgo'])
print(table.round(3))
ecuación del libro descenso del gradiente mínimos cuadrados \
pacientes 0.05 0.044 0.044
temperatura -0.37 -0.374 -0.374
urgencias 3.20 3.207 3.207
alerta 18.70 18.808 18.808
fin_de_semana 4.10 3.363 3.363
sesgo 42.30 44.304 44.304
scikit-learn
pacientes 0.044
temperatura -0.374
urgencias 3.207
alerta 18.808
fin_de_semana 3.363
sesgo 44.304
Los tres métodos coinciden hasta el tercer decimal y se acercan a la ecuación del libro: 0,044 frente a 0,05 por paciente, −0,374 frente a −0,37 por grado, 3,207 frente a 3,2 por urgencia, 18,808 frente a 18,7 por alerta y 44,304 frente a 42,3 de sesgo. La mayor diferencia es la del fin de semana, 3,363 frente a 4,1: solo uno de cada siete días es de fin de semana, así que su peso se estima con menos días y el ruido se nota más en él. Con datos reales no hay ecuación con la que comparar, y por eso el paso siguiente prueba el modelo con días que no ha visto.
Paso 8: ¿Cómo de bueno es el modelo?¶
Ahora usamos los tres meses de prueba, que ningún método ha visto. Además del ECM y del EAM (el error absoluto medio, en mascarillas), calculamos el coeficiente de determinación $R^2 = 1 - \frac{\sum (y - \hat{y})^2}{\sum (y - \bar{y})^2}$, la parte de la variación de las mascarillas de un día a otro que explica el modelo.
y_pred = model.predict(X_test)
print(f"ECM de prueba: {mean_squared_error(y_test, y_pred):.2f} (el ruido que añadimos tiene varianza 4² = 16)")
print(f"EAM de prueba: {mean_absolute_error(y_test, y_pred):.2f} mascarillas")
print(f"R² de prueba: {r2_score(y_test, y_pred):.3f}")
plt.figure(figsize=(12, 4))
plt.plot(dates[~train], y_test, label='mascarillas usadas')
plt.plot(dates[~train], y_pred, label='predicción')
plt.title('Los tres meses de prueba')
plt.legend()
plt.show()
ECM de prueba: 13.34 (el ruido que añadimos tiene varianza 4² = 16) EAM de prueba: 2.75 mascarillas R² de prueba: 0.889
En los tres meses de prueba, el ECM es 13,34, el EAM 2,75 mascarillas y el R² 0,889: el modelo explica cerca del 89% de la variación de un día a otro, y lo que queda es sobre todo el ruido que pusimos en los datos. La predicción sigue de cerca las mascarillas usadas.
Paso 9: Cuando los datos tienen errores: ¿ECM o EAM?¶
En el libro, el equipo compara cuatro funciones de coste y elige el ECM, con una regla de Luis: si aparecen valores atípicos, pasar al EAM. El error al cuadrado hace que un fallo grande cuente mucho más (un error de 10 cuesta 100, como señala Marta), así que unos pocos registros erróneos pueden arrastrar todo el modelo. Añadimos cinco errores de tecleo a los días de entrenamiento, mascarillas escritas con un cero de más, y ajustamos el modelo dos veces: minimizando el ECM (LinearRegression) y minimizando el EAM (QuantileRegressor con la mediana, quantile=0.5, y sin penalización).
y_typos = y_train.copy()
wrong_days = [40, 130, 222, 300, 390]
y_typos[wrong_days] = y_train[wrong_days] * 10 # p. ej., 85 escrito como 850
print("Cinco registros erróneos:", y_typos[wrong_days].round())
fits = {'ECM, datos limpios': LinearRegression().fit(X_train, y_train),
'ECM, con los errores': LinearRegression().fit(X_train, y_typos),
'EAM, con los errores': QuantileRegressor(quantile=0.5, alpha=0, solver='highs').fit(X_train, y_typos)}
rows = {name: np.r_[f.coef_, f.intercept_, mean_absolute_error(y_test, f.predict(X_test))] for name, f in fits.items()}
print(pd.DataFrame(rows, index=features + ['sesgo', 'EAM de prueba']).round(3))
Cinco registros erróneos: [1146. 934. 966. 813. 1106.]
ECM, datos limpios ECM, con los errores EAM, con los errores
pacientes 0.044 0.220 0.039
temperatura -0.374 -0.504 -0.390
urgencias 3.207 4.601 3.167
alerta 18.808 4.373 18.741
fin_de_semana 3.363 36.153 3.234
sesgo 44.304 -17.253 46.452
EAM de prueba 2.750 11.607 2.755
Cinco registros erróneos de 456 bastan para estropear el ajuste con ECM: el peso de la alerta baja de 18,808 a 4,373, el del fin de semana se dispara a 36,153 y el EAM de prueba sube de 2,75 a 11,607 mascarillas. El ajuste con EAM, entrenado con los mismos datos erróneos, se queda casi donde estaba el ajuste limpio (alerta 18,741, EAM de prueba 2,755). Es la regla de Luis en números. Aun así, el primer remedio es encontrar y corregir los registros erróneos: el EAM limita el daño, no nos dice dónde están los errores.
Paso 10: La previsión para la semana que viene¶
Por último, la semana con la que el equipo prueba su modelo en el libro: el lunes, con 295 pacientes, 14 °C y 12 urgencias previstas; el martes, con 310 pacientes, 15 °C y 10 urgencias; el miércoles, con 285 pacientes, 16 °C y 11 urgencias; sin alerta sanitaria y sin fin de semana.
week = pd.DataFrame([[295, 14, 12, 0, 0], [310, 15, 10, 0, 0], [285, 16, 11, 0, 0]],
index=['Lunes', 'Martes', 'Miércoles'], columns=features)
week['ecuación del libro'] = week[features].values @ book_w + book_b
week['nuestro modelo'] = model.predict(week[features].values)
print(week.round(2))
pacientes temperatura urgencias alerta fin_de_semana \
Lunes 295 14 12 0 0
Martes 310 15 10 0 0
Miércoles 285 16 11 0 0
ecuación del libro nuestro modelo
Lunes 90.27 90.63
Martes 84.25 84.51
Miércoles 85.83 86.23
Nuestro modelo predice 90,63, 84,51 y 86,23 mascarillas para el lunes, el martes y el miércoles. La ecuación del libro da exactamente 90,27, 84,25 y 85,83, que el libro redondea hacia arriba a 91, 85 y 86. Las dos previsiones están dentro de los rangos históricos de Carlos: de 85 a 95 mascarillas para un lunes como este, de 82 a 86 para el martes y de 80 a 90 para el miércoles.
Conclusión¶
- Con varias variables, la regresión lineal sigue siendo una suma ponderada, $\hat{y} = \mathbf{w} \cdot \mathbf{x} + b$: un hiperplano en lugar de una recta, y cada peso se lee como el efecto de una variable cuando las demás no cambian.
- Las ecuaciones de actualización de Luis funcionan igual con cualquier número de variables, pero solo si las variables están en escalas parecidas. Estandarizarlas convierte miles de iteraciones inútiles en unos cientos útiles.
- Para la regresión lineal, los mínimos cuadrados dan la respuesta exacta directamente, y scikit-learn los usa. El descenso del gradiente importa porque es el método que funciona también con los modelos de los capítulos siguientes.
- La función de coste es una decisión: el ECM es la elección habitual, pero unos pocos registros erróneos pueden estropearlo, y el EAM los resiste.