Análisis Avanzado de Datos — 5. Validación cruzada k-fold

estadistica
validacion-cruzada
k-fold
Complemento de evaluación: validación cruzada k-fold paso a paso, particiones, elección de k e interpretación en regresión lineal, con gráficos de diagnóstico.
Author

Wilder Ramírez Delgado

Published

September 8, 2026

Análisis Avanzado de Datos — 5. Validación cruzada k-fold

Open in Colab

👋 Sobre el autor

Wilder Ramírez Delgado convierte datos en decisiones. Científico de Datos, Arquitecto de IA y asesor tecnológico con más de 20 años de experiencia en analítica avanzada, machine learning, Big Data e IoT. Ingeniero Electrónico y Magíster en Analítica de Datos. CEO y fundador de Business Innovation Technology (BIT), asesora a entidades gubernamentales y empresas privadas en Colombia y es referente en analítica y normatividad para los sectores postal y logístico del país. Una serie pensada para ir de la teoría a la práctica, con ejemplos concretos en cada paso.

📌 Parte de la serie Análisis Avanzado de Datos.

Un solo número de prueba es una lotería

Entrenaste un modelo, lo evaluaste en un conjunto de prueba y obtuviste un \(R^2\) de 0.61. Lo presentas. Un colega toma los mismos datos, hace el mismo modelo y obtiene 0.39. ¿Quién midió mal?

Ninguno. Cada uno reservó un grupo distinto de observaciones para la prueba. El desempeño que reportas no es una propiedad fija del modelo: es una realización que depende de qué filas cayeron de un lado o del otro. Con una sola partición estás apostando a una carta.

La validación cruzada k-fold (validación cruzada en \(k\) bloques) es la forma sistemática de dejar de apostar. En vez de una sola división, rotas qué bloque queda fuera, mides \(k\) veces y reportas un promedio acompañado de su dispersión. El objetivo de este notebook no es solo ejecutar cross_val_score: es que veas las particiones, entiendas el proceso observación por observación, y después sepas leer los números que salen cuando el modelo es una recta y cuando el modelo tiene varias variables.

Al terminar vas a poder:

  • Explicar, con un dibujo, qué hace k-fold y por qué cada observación se usa exactamente una vez como prueba.
  • Distinguir lo que k-fold estima (el error de predicción esperado del procedimiento) de lo que no estima (el error de un modelo ya congelado).
  • Elegir \(k\) con criterio y detectar el error de no barajar cuando los datos vienen ordenados.
  • Leer un resultado fold a fold: no solo el promedio, sino por qué un bloque salió peor.
  • Aplicarlo a una regresión lineal simple y a una regresión multivariable, con gráficos de diagnóstico en cada paso.

Preparación del entorno

NumPy, pandas, matplotlib, scikit-learn y statsmodels, más una semilla fija (SEMILLA) para que los gráficos y los números salgan iguales cada vez que corras el notebook. Los datos de esta primera parte son sintéticos y pequeños a propósito: así se puede dibujar cada observación.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.colors import ListedColormap

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.datasets import load_diabetes
import statsmodels.api as sm

SEMILLA = 2026
rng = np.random.default_rng(SEMILLA)

plt.rcParams['figure.dpi'] = 110
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False

COLOR_TRAIN = '#4C78A8'
COLOR_TEST = '#E45756'
COLOR_LINEA = '#333333'
colores = ['#4C78A8', '#F58518', '#E45756', '#54A24B', '#B279A2']

El problema que k-fold viene a resolver

Toma un dataset sintético pequeño, pártelo 70/30 de 40 formas distintas y guarda el \(R^2\) de prueba de cada partición. Si el número “del modelo” fuera estable, esos 40 valores se parecerían y el histograma se vería como un clavo. No lo es.

n_lot = 80
x_lot = rng.uniform(0, 10, size=n_lot)
y_lot = 2.0 + 1.5 * x_lot + rng.normal(0, 2.2, size=n_lot)

r2_splits = []
for s in range(40):
    Xa, Xb, ya, yb = train_test_split(
        x_lot.reshape(-1, 1), y_lot, test_size=0.3, random_state=s
    )
    m = LinearRegression().fit(Xa, ya)
    r2_splits.append(r2_score(yb, m.predict(Xb)))

r2_splits = np.array(r2_splits)

print(f'R2 mínimo:   {r2_splits.min():.3f}')
print(f'R2 máximo:   {r2_splits.max():.3f}')
print(f'R2 promedio: {r2_splits.mean():.3f}')
print(f'Desviación:  {r2_splits.std():.3f}')
print(f'Rango:       {r2_splits.max() - r2_splits.min():.3f}')
R2 mínimo:   0.577
R2 máximo:   0.831
R2 promedio: 0.725
Desviación:  0.065
Rango:       0.254

Ver la nube de puntos y el histograma de esos 40 R²

A la izquierda, los datos. A la derecha, la lotería: cómo cambia el \(R^2\) de prueba según qué filas cayeron fuera.

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

ax = axes[0]
ax.scatter(x_lot, y_lot, c=COLOR_TRAIN, s=28, alpha=0.85, edgecolors='white', linewidths=0.4)
xs = np.linspace(0, 10, 50)
ax.plot(xs, 2.0 + 1.5 * xs, color=COLOR_LINEA, lw=1.5, label='Relación verdadera')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Los datos (n = 80)')
ax.legend(frameon=False)

ax = axes[1]
ax.hist(r2_splits, bins=12, color=COLOR_TRAIN, edgecolor='white')
ax.axvline(r2_splits.mean(), color=COLOR_TEST, ls='--', lw=1.5,
           label=f'Promedio = {r2_splits.mean():.3f}')
ax.set_xlabel('R² de prueba (una partición 70/30)')
ax.set_ylabel('Frecuencia')
ax.set_title('40 formas distintas de cortar los mismos datos')
ax.legend(frameon=False)

plt.tight_layout()
plt.show()

Cómo leer este resultado

A la izquierda no hay ningún misterio: una nube lineal con ruido. A la derecha, el mismo modelo y los mismos 80 puntos producen un \(R^2\) de prueba que va de 0.58 a 0.83 (rango de 0.25) según qué 24 observaciones cayeron en la prueba.

Si te hubiera tocado la partición de 0.83, reportarías un modelo claramente mejor de lo que es. Si te tocara el de 0.58, lo subestimarías. Ninguno de los dos extremos es “el” desempeño — el centro está cerca de 0.73. k-fold existe para estimar ese centro de forma ordenada —y para mostrarte, fold a fold, cuánto varía el error—, no con 40 particiones improvisadas.

Qué es k-fold, en una frase y en fórmulas

Divides las \(n\) observaciones en \(k\) bloques (folds) de tamaño parecido. En la ronda \(j\) entrenas con todos los bloques menos el \(j\) y evalúas solo en el \(j\). Rotas \(j = 1, \dots, k\). Al final, cada observación fue prueba exactamente una vez y entrenamiento en las otras \(k-1\) rondas.

Si \(e_j\) es el error (o el \(R^2\)) medido en el fold \(j\), el estimador de validación cruzada es el promedio:

\[\mathrm{CV}_k = \frac{1}{k}\sum_{j=1}^{k} e_j\]

y la dispersión entre folds se resume con

\[s_{\mathrm{CV}} = \sqrt{\frac{1}{k}\sum_{j=1}^{k}\bigl(e_j - \mathrm{CV}_k\bigr)^2}\]

Dos precisiones que importan más de lo que parecen:

  1. \(\mathrm{CV}_k\) estima el error de predicción esperado del procedimiento (el algoritmo + el tamaño de muestra de entrenamiento), no el error futuro de un modelo ya ajustado con el 100% de los datos. Son primos, no gemelos.
  2. \(s_{\mathrm{CV}}\) describe cuánto varía el error entre bloques. No es, en sentido estricto, el error estándar de \(\mathrm{CV}_k\): los folds no son independientes, porque comparten observaciones de entrenamiento. Úsalo como termómetro de estabilidad, no como un error estándar para construir intervalos de confianza.

El valor por defecto en la práctica es \(k=5\) o \(k=10\). Con \(k=5\) cada modelo se entrena con el 80% de los datos; con \(k=10\), con el 90%. \(k=n\) es el caso extremo (leave-one-out): entrenas \(n\) veces dejando una sola observación fuera. Casi no tiene sesgo por tamaño de entrenamiento, pero su varianza puede ser mayor y exige \(n\) ajustes (en regresión lineal existe un atajo exacto, pero no en general).

Las particiones, dibujadas

Antes del código de scikit-learn, mira quién queda fuera en cada ronda. Con \(n=20\) y \(k=5\) cada fold de prueba tiene 4 observaciones: lo suficientemente pequeño para que puedas seguir cada punto. El código imprime, fold a fold, el índice de las observaciones de prueba y comprueba que cada una aparece exactamente una vez.

n_demo = 20
k = 5
x_demo = np.sort(rng.uniform(1, 10, size=n_demo))
y_demo = 2.0 + 1.5 * x_demo + rng.normal(0, 1.1, size=n_demo)

kf_demo = KFold(n_splits=k, shuffle=True, random_state=SEMILLA)

matriz = np.zeros((k, n_demo))
for j, (_, te) in enumerate(kf_demo.split(x_demo)):
    matriz[j, te] = 1
    print(f'Fold {j+1}: prueba en observaciones {[int(i) + 1 for i in te]} '
          f'(n_train={n_demo - len(te)}, n_test={len(te)})')

print()
print('Comprobación: cada observación aparece exactamente una vez como prueba:')
print('veces en prueba =', matriz.sum(axis=0).astype(int))
Fold 1: prueba en observaciones [10, 11, 12, 16] (n_train=16, n_test=4)
Fold 2: prueba en observaciones [8, 9, 15, 20] (n_train=16, n_test=4)
Fold 3: prueba en observaciones [1, 3, 6, 19] (n_train=16, n_test=4)
Fold 4: prueba en observaciones [4, 5, 13, 17] (n_train=16, n_test=4)
Fold 5: prueba en observaciones [2, 7, 14, 18] (n_train=16, n_test=4)

Comprobación: cada observación aparece exactamente una vez como prueba:
veces en prueba = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

Dibujar la cuadrícula de particiones

Cada fila es una ronda; azul es entrenamiento y rojo es prueba. Si el contrato de k-fold se cumple, en cada columna hay un solo rojo.

fig, ax = plt.subplots(figsize=(11, 3.4))
cmap = ListedColormap([COLOR_TRAIN, COLOR_TEST])
ax.imshow(matriz, aspect='auto', cmap=cmap, vmin=0, vmax=1, interpolation='nearest')
ax.set_xticks(range(n_demo))
ax.set_xticklabels([str(i) for i in range(1, n_demo + 1)])
ax.set_yticks(range(k))
ax.set_yticklabels([f'Fold {j+1}' for j in range(k)])
ax.set_xlabel('Observación (índice)')
ax.set_title('Particiones k-fold (k = 5, n = 20): azul = entrena, rojo = prueba')
ax.legend(
    handles=[
        Patch(facecolor=COLOR_TRAIN, label='Entrenamiento'),
        Patch(facecolor=COLOR_TEST, label='Prueba'),
    ],
    loc='upper left',
    bbox_to_anchor=(1.01, 1.0),
    frameon=False,
)
plt.tight_layout()
plt.show()

Los cuadrados rojos nunca se superponen en la misma columna: si la observación 7 es prueba en el fold 5, en los otros cuatro folds es entrenamiento. Eso es el contrato de k-fold.

Ver el proceso: ajustar con el azul y medir el residual del rojo

La siguiente figura es el proceso, no solo la etiqueta: cinco paneles, uno por fold. En cada uno la recta se ajusta solo con los puntos azules y el error se mide solo en los rojos. Las rayas verticales son los residuales de prueba: la distancia entre lo que predice el modelo, que no vio ese punto, y lo que realmente se midió.

fig, axes = plt.subplots(1, k, figsize=(13.5, 3.3), sharex=True, sharey=True)
xs = np.linspace(x_demo.min() - 0.3, x_demo.max() + 0.3, 60)

for j, (tr, te) in enumerate(kf_demo.split(x_demo)):
    ax = axes[j]
    m = LinearRegression().fit(x_demo[tr].reshape(-1, 1), y_demo[tr])
    pred_te = m.predict(x_demo[te].reshape(-1, 1))
    ax.scatter(x_demo[tr], y_demo[tr], c=COLOR_TRAIN, s=32, zorder=2, label='Entrena')
    ax.scatter(x_demo[te], y_demo[te], c=COLOR_TEST, s=55, zorder=3, label='Prueba')
    ax.plot(xs, m.predict(xs.reshape(-1, 1)), color=COLOR_LINEA, lw=1.3)
    for xj, yj, pj in zip(x_demo[te], y_demo[te], pred_te):
        ax.plot([xj, xj], [pj, yj], color=COLOR_TEST, lw=1.0, alpha=0.85)
    rmse = np.sqrt(mean_squared_error(y_demo[te], pred_te))
    ax.set_title(f'Fold {j+1}\nRMSE = {rmse:.2f}')
    ax.set_xlabel('x')
    if j == 0:
        ax.set_ylabel('y')
        ax.legend(frameon=False, fontsize=8, loc='upper left')

fig.suptitle('El proceso: ajustar con el azul, evaluar el residual del rojo', y=1.05)
plt.tight_layout()
plt.show()

Fíjate en tres cosas a la vez:

  • La recta cambia de un fold a otro, porque cada ronda deja fuera 4 puntos distintos. Eso es normal: estás viendo la variabilidad del estimador.
  • El RMSE de prueba (raíz del error cuadrático medio, en las unidades de \(y\); lo vimos en el artículo de evaluación) no es idéntico en cada panel: va de 0.81 (fold 2) a 2.34 (fold 5). Un fold sale peor cuando sus puntos rojos son ruidosos: en el fold 5, la observación con \(x \approx 7.5\) queda más de 4 unidades por debajo de su recta; en el fold 3, la de \(x \approx 1.3\), unas 3 unidades.
  • Ningún punto rojo se usó para trazar su propia recta. Por eso este error es más honesto que el error de entrenamiento.

Esa rotación —cinco modelos, cinco pruebas, un promedio— es todo k-fold. El resto de este notebook es aprender a leerla cuando los datos dejan de caber en veinte puntos.

Por qué hay que barajar (salvo que el orden importe)

KFold por defecto no baraja: parte por orden de índice, en bloques consecutivos. Si tus filas están ordenadas por tiempo, por tamaño o por el valor de \(y\), cada fold deja de ser una muestra del mismo fenómeno. Uno se queda con los valores pequeños de \(x\); otro, con los grandes. Entonces no estás midiendo estabilidad del modelo: estás midiendo que los extremos no se parecen al centro.

Cuando las observaciones se pueden considerar intercambiables (el caso habitual en corte transversal), usa shuffle=True y fija random_state para que el resultado sea reproducible. Cuando el orden es la información —series de tiempo, mediciones repetidas del mismo paciente— k-fold aleatorio deja de ser válido y hace falta otro diseño (bloques temporales, agrupado). Aquí nos quedamos en el caso i.i.d. (observaciones independientes e idénticamente distribuidas).

Comparar k-fold sin barajar y con shuffle

Los datos van ordenados de menor a mayor \(x\). Sin shuffle, cada fold de prueba es una franja del eje; con shuffle, cada color se mezcla a lo largo de todo el rango.

n_ord = 40
x_ord = np.linspace(1, 10, n_ord)
y_ord = 2.0 + 1.5 * x_ord + rng.normal(0, 1.0, size=n_ord)

kf_no = KFold(n_splits=5, shuffle=False)
kf_si = KFold(n_splits=5, shuffle=True, random_state=SEMILLA)

fold_no = np.empty(n_ord, dtype=int)
fold_si = np.empty(n_ord, dtype=int)
for j, (_, te) in enumerate(kf_no.split(x_ord)):
    fold_no[te] = j + 1
for j, (_, te) in enumerate(kf_si.split(x_ord)):
    fold_si[te] = j + 1

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), sharex=True, sharey=True)
for ax, folds, titulo in [
    (axes[0], fold_no, 'Sin barajar: los folds son franjas de x'),
    (axes[1], fold_si, 'Con shuffle: cada fold mezcla todo el rango'),
]:
    for j in range(1, 6):
        mask = folds == j
        ax.scatter(x_ord[mask], y_ord[mask], c=colores[j - 1], s=36,
                   label=f'Fold {j}', edgecolors='white', linewidths=0.4)
    ax.set_xlabel('x (ordenado de menor a mayor)')
    ax.set_title(titulo)
    ax.legend(frameon=False, fontsize=8, loc='upper left')
axes[0].set_ylabel('y')
plt.tight_layout()
plt.show()

print('Media de x en cada fold de prueba')
resumen_fold = pd.DataFrame({
    'fold': np.arange(1, 6),
    'media_x_sin_barajar': [x_ord[fold_no == j].mean() for j in range(1, 6)],
    'media_x_con_shuffle': [x_ord[fold_si == j].mean() for j in range(1, 6)],
})
print(resumen_fold.round(2).to_string(index=False))

Media de x en cada fold de prueba
 fold  media_x_sin_barajar  media_x_con_shuffle
    1                 1.81                 6.28
    2                 3.65                 4.58
    3                 5.50                 6.57
    4                 7.35                 4.92
    5                 9.19                 5.15

A la izquierda, el fold 1 prueba solo la zona baja de \(x\) y el fold 5 solo la zona alta. En esos dos folds le estás pidiendo al modelo que extrapole, así que su error mide otra cosa que el de los folds centrales. Aquí la relación verdadera es lineal y el costo es pequeño; con un modelo mal especificado (una curva ajustada con una recta), extrapolar en los extremos dispara el error. A la derecha, cada color está mezclado a lo largo de todo el eje. Esa es la partición que quieres cuando asumes que las filas son intercambiables.

La tabla lo dice en números: sin barajar, la media de \(x\) en prueba camina de forma monótona; con shuffle, las cinco medias se parecen.

Cómo elegir \(k\) (y qué no esperar de esa elección)

No hay un \(k\) mágico. El compromiso es este:

  • \(k\) pequeño (2 o 3): entrenas con pocos datos, el modelo de cada ronda es peor que el que ajustarías con \(n\), y \(\mathrm{CV}_k\) tiende a ser un poco pesimista. Barato de computar.
  • \(k=5\) o \(k=10\): el estándar. Entrenamiento al 80–90%, sesgo pequeño, varianza razonable, costo aceptable.
  • \(k=n\) (leave-one-out): cada modelo se entrena con \(n-1\) filas, así que el sesgo de tamaño casi desaparece. A cambio, las \(n\) medidas están muy correlacionadas y el estimador puede tener más varianza; además entrenas \(n\) veces.

Para \(n\) de unos cientos, \(k=5\) es una elección seria, no una pereza. Más adelante, en los ejemplos, vas a comparar \(k=5\) y \(k=10\) sobre los mismos datos para ver cuánto se mueve el número — adelanto: mucho menos de lo que se movía una sola partición 70/30 en la primera sección.

Un error frecuente: usar k-fold para elegir el modelo (o un hiperparámetro) y reportar ese mismo \(\mathrm{CV}_k\) como desempeño final. Ese número queda optimista. La versión honesta es validación cruzada anidada, o una prueba que no haya participado en la elección. En este notebook el modelo está fijo —una regresión lineal—, así que un solo k-fold sí es una estimación honesta de su error de predicción.

Ejemplo 1. Regresión lineal simple: precio a partir de los metros

Hasta aquí el dibujo. Ahora un problema con unidades. Imagina que tasas apartamentos y, en una primera pasada, solo tienes los metros cuadrados para predecir el precio (en miles de dólares). La relación verdadera —que en la vida real no conoces, y aquí sí, porque los datos son simulados— es

\[\mathrm{precio} = 28 + 1.80 \cdot \mathrm{metros} + \varepsilon\]

con un ruido que crece un poco en los apartamentos grandes y con tres tasaciones atípicas. k-fold tiene que convivir con eso: no con una nube de laboratorio. El código simula 120 apartamentos e imprime un describe, los tres atípicos y la correlación entre metros y precio.

n1 = 120
metros = rng.uniform(45, 165, size=n1)
sigma = 9 + 0.07 * (metros - 45)
precio = 28 + 1.80 * metros + rng.normal(0, 1, size=n1) * sigma

idx_atipicos = np.array([18, 61, 104])
precio[idx_atipicos] = precio[idx_atipicos] + np.array([52.0, -48.0, 58.0])

df1 = pd.DataFrame({
    'metros': metros,
    'precio': precio,
    'atipico': np.isin(np.arange(n1), idx_atipicos),
})

print(df1[['metros', 'precio']].describe().round(2).to_string())
print()
print('Apartamentos atípicos (índice, metros, precio):')
print(df1.loc[idx_atipicos, ['metros', 'precio']].round(1).to_string())
print()
print(f'Correlación metros–precio: {df1["metros"].corr(df1["precio"]):.3f}')
       metros  precio
count  120.00  120.00
mean   109.14  225.59
std     35.32   65.65
min     46.11  105.50
25%     82.67  173.71
50%    107.87  223.10
75%    138.92  273.01
max    164.93  354.44

Apartamentos atípicos (índice, metros, precio):
     metros  precio
18    150.1   336.4
61    146.0   237.2
104    47.9   166.1

Correlación metros–precio: 0.969

Explorar metros, precio y la nube con la recta verdadera

Tres gráficos: histograma del predictor, histograma del objetivo, y scatter con la recta que sí conocemos (porque los datos son simulados) y los atípicos marcados.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.0))

ax = axes[0]
ax.hist(df1['metros'], bins=15, color=COLOR_TRAIN, edgecolor='white')
ax.set_xlabel('Metros cuadrados')
ax.set_ylabel('Frecuencia')
ax.set_title('Predictor: metros')

ax = axes[1]
ax.hist(df1['precio'], bins=15, color='#54A24B', edgecolor='white')
ax.set_xlabel('Precio (miles de USD)')
ax.set_title('Objetivo: precio')

ax = axes[2]
normales = ~df1['atipico']
ax.scatter(df1.loc[normales, 'metros'], df1.loc[normales, 'precio'],
           c=COLOR_TRAIN, s=28, alpha=0.85, edgecolors='white', linewidths=0.4,
           label='Observaciones')
ax.scatter(df1.loc[df1['atipico'], 'metros'], df1.loc[df1['atipico'], 'precio'],
           c=COLOR_TEST, s=70, marker='D', zorder=3, label='Atípicos')
xs = np.linspace(45, 165, 50)
ax.plot(xs, 28 + 1.80 * xs, color=COLOR_LINEA, lw=1.5, label='Relación verdadera')
ax.set_xlabel('Metros cuadrados')
ax.set_ylabel('Precio (miles de USD)')
ax.set_title('Nube, recta verdadera y 3 atípicos')
ax.legend(frameon=False, fontsize=8)

plt.tight_layout()
plt.show()

Cómo leer estos datos

Los metros cubren de ~45 a ~165 m², sin huecos raros. El precio se reparte de forma bastante pareja entre ~105 y ~355 mil USD, sin una cola marcada. Los atípicos que marcamos en rojo no son extremos del precio (336, 237 y 166 mil USD: uno arriba, otro en la mitad y otro en la parte baja del rango): son atípicos dado su tamaño, es decir, se alejan de la recta, no del histograma. En el diagrama de dispersión la tendencia es lineal y clara —la correlación sale 0.97—, pero el abanico se abre a la derecha: los apartamentos grandes tienen más ruido. Los tres diamantes no son un error de código; son tasaciones que un analista vería como “esta se salió”. k-fold los va a poner, en alguna ronda, del lado de la prueba. Eso es bueno: quieres saber cuánto duele un punto así cuando el modelo no lo ha visto.

Agregar el intercepto y ajustar OLS con toda la muestra

Antes de rotar folds, ajusta la recta con todos los datos. X1 son solo los metros. statsmodels no pone ordenada al origen solo: add_constant agrega una columna de unos, y con esa matriz se estima la recta por mínimos cuadrados ordinarios (OLS, por sus siglas en inglés). El summary() describe el ajuste in-sample (sobre los mismos datos con que se estimó), todavía no k-fold: útil como referencia, optimista como evaluación.

X1 = df1[['metros']].values
y1 = df1['precio'].values

X1_sm = sm.add_constant(X1)
ols1 = sm.OLS(y1, X1_sm).fit()
print(ols1.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.939
Model:                            OLS   Adj. R-squared:                  0.938
Method:                 Least Squares   F-statistic:                     1811.
Date:                Sun, 27 Sep 2026   Prob (F-statistic):           1.90e-73
Time:                        20:55:50   Log-Likelihood:                -504.24
No. Observations:                 120   AIC:                             1012.
Df Residuals:                     118   BIC:                             1018.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         29.0075      4.853      5.978      0.000      19.398      38.617
x1             1.8011      0.042     42.562      0.000       1.717       1.885
==============================================================================
Omnibus:                        3.977   Durbin-Watson:                   2.032
Prob(Omnibus):                  0.137   Jarque-Bera (JB):                4.730
Skew:                           0.086   Prob(JB):                       0.0940
Kurtosis:                       3.957   Cond. No.                         374.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

El intercepto y la pendiente deberían acercarse a 28 y 1.80. No van a caer encima: hay ruido y tres atípicos. El \(R^2\) de esta salida es in-sample: dice qué tanto la recta describe estos 120 puntos, no qué tanto predeciría 120 apartamentos nuevos. Eso lo mide k-fold.

Diagnosticar ese ajuste in-sample

Tres scatters: la recta estimada contra la verdadera, predicho vs observado, y residuales vs predicho. Sirven de referencia; el veredicto honesto viene con los folds.

pred_in = ols1.fittedvalues
resid_in = ols1.resid

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.0))

ax = axes[0]
ax.scatter(df1['metros'], df1['precio'], c=COLOR_TRAIN, s=28, alpha=0.8,
           edgecolors='white', linewidths=0.4)
xs = np.linspace(df1['metros'].min(), df1['metros'].max(), 50)
ax.plot(xs, ols1.params[0] + ols1.params[1] * xs, color=COLOR_TEST, lw=1.6,
        label='OLS (muestra completa)')
ax.plot(xs, 28 + 1.80 * xs, color=COLOR_LINEA, lw=1.2, ls='--',
        label='Verdadera')
ax.set_xlabel('Metros')
ax.set_ylabel('Precio')
ax.set_title('Ajuste in-sample')
ax.legend(frameon=False, fontsize=8)

ax = axes[1]
ax.scatter(pred_in, y1, c=COLOR_TRAIN, s=28, alpha=0.8,
           edgecolors='white', linewidths=0.4)
lo, hi = pred_in.min(), pred_in.max()
ax.plot([lo, hi], [lo, hi], color=COLOR_LINEA, lw=1.2)
ax.set_xlabel('Precio predicho')
ax.set_ylabel('Precio observado')
ax.set_title('Predicho vs observado (entrenamiento)')

ax = axes[2]
ax.scatter(pred_in, resid_in, c=COLOR_TRAIN, s=28, alpha=0.8,
           edgecolors='white', linewidths=0.4)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('Precio predicho')
ax.set_ylabel('Residual')
ax.set_title('Residuales in-sample')

plt.tight_layout()
plt.show()

print(f'Pendiente OLS: {ols1.params[1]:.3f}  (verdadera = 1.80)')
print(f'Intercepto OLS: {ols1.params[0]:.2f}  (verdadero = 28)')
print(f'R2 in-sample: {ols1.rsquared:.3f}')
print(f'RMSE in-sample: {np.sqrt(np.mean(resid_in**2)):.2f}')

Pendiente OLS: 1.801  (verdadera = 1.80)
Intercepto OLS: 29.01  (verdadero = 28)
R2 in-sample: 0.939
RMSE in-sample: 16.17

La recta estimada pasa cerca de la verdadera. El gráfico predicho vs observado se pega a la diagonal, con tres puntos que se salen. Los residuales se abren un poco a la derecha: heterocedasticidad suave, la que inyectamos. Nada de esto es todavía una evaluación honesta. Ahora sí: k-fold manual, para ver cada bloque.

Recorrer los 5 folds a mano y tabular el resultado de cada bloque

En cada ronda se entrena con el 80% y se evalúa en el 20% que sobra. La tabla guarda, por fold, rango de metros, intercepto, pendiente, MAE (error absoluto medio), RMSE, \(R^2\) y cuántos atípicos cayeron en prueba.

k5 = KFold(n_splits=5, shuffle=True, random_state=SEMILLA)
fold_id = np.empty(n1, dtype=int)

filas = []
pred_cv = np.empty(n1)
resid_cv = np.empty(n1)

for j, (tr, te) in enumerate(k5.split(X1)):
    fold_id[te] = j + 1
    m = LinearRegression().fit(X1[tr], y1[tr])
    pred = m.predict(X1[te])
    pred_cv[te] = pred
    resid_cv[te] = y1[te] - pred
    filas.append({
        'fold': j + 1,
        'n_test': len(te),
        'metros_promedio': X1[te, 0].mean(),
        'metros_min': X1[te, 0].min(),
        'metros_max': X1[te, 0].max(),
        'precio_promedio': y1[te].mean(),
        'precio_sd': y1[te].std(),
        'intercepto': m.intercept_,
        'pendiente': m.coef_[0],
        'MAE': mean_absolute_error(y1[te], pred),
        'RMSE': np.sqrt(mean_squared_error(y1[te], pred)),
        'R2': r2_score(y1[te], pred),
        'n_atipicos_en_prueba': int(np.isin(te, idx_atipicos).sum()),
    })

tab1 = pd.DataFrame(filas)
print(tab1.round(2).to_string(index=False))
print()
print(f'RMSE promedio (CV): {tab1["RMSE"].mean():.2f}  ± {tab1["RMSE"].std(ddof=0):.2f}')
print(f'R2 promedio   (CV): {tab1["R2"].mean():.3f}  ± {tab1["R2"].std(ddof=0):.3f}')
print(f'Pendiente media:    {tab1["pendiente"].mean():.3f}  '
      f'(rango {tab1["pendiente"].min():.3f} a {tab1["pendiente"].max():.3f})')
 fold  n_test  metros_promedio  metros_min  metros_max  precio_promedio  precio_sd  intercepto  pendiente   MAE  RMSE   R2  n_atipicos_en_prueba
    1      24           118.34       46.11      164.27           242.50      64.94       23.07       1.86 17.00 22.65 0.88                     3
    2      24           112.22       62.65      162.13           228.58      53.59       30.60       1.79 12.39 16.04 0.91                     0
    3      24           108.55       51.06      164.93           224.00      69.97       33.35       1.76 10.88 14.14 0.96                     0
    4      24           115.20       46.98      163.03           238.31      70.66       32.08       1.77 12.29 14.13 0.96                     0
    5      24            91.41       46.23      158.06           194.54      54.50       25.38       1.83 11.56 14.68 0.93                     0

RMSE promedio (CV): 16.33  ± 3.24
R2 promedio   (CV): 0.927  ± 0.031
Pendiente media:    1.802  (rango 1.763 a 1.856)

Pintar cada apartamento con el fold en el que fue prueba

Los atípicos van con borde negro para ver en qué bloque cayeron.

fig, ax = plt.subplots(figsize=(8, 4.5))
for j in range(1, 6):
    mask = fold_id == j
    ax.scatter(df1.loc[mask, 'metros'], df1.loc[mask, 'precio'],
               c=colores[j - 1], s=32, label=f'Fold {j}',
               edgecolors='white', linewidths=0.4, zorder=2)
ax.scatter(df1.loc[idx_atipicos, 'metros'], df1.loc[idx_atipicos, 'precio'],
           facecolors='none', edgecolors='black', s=90, lw=1.2, zorder=3,
           label='Atípicos')
ax.set_xlabel('Metros cuadrados')
ax.set_ylabel('Precio (miles de USD)')
ax.set_title('Qué apartamentos cayeron en cada fold de prueba')
ax.legend(frameon=False, fontsize=8, loc='upper left')
plt.tight_layout()
plt.show()

El scatter coloreado te dice si el shuffle funcionó: cada color debe cubrir el rango de metros, no un pedazo. Los círculos negros son los atípicos; anota en qué fold de prueba cayeron (columna n_atipicos_en_prueba de la tabla). El fold que se los lleve debería mostrar un RMSE más alto o un \(R^2\) más bajo. Si eso ocurre, ya no es “ruido misterioso del CV”: es un mecanismo.

Ver las cinco rectas: cada una ignora su bloque rojo

Mismo scatter, cinco veces. Los rojos no participaron en trazar esa recta. La pendiente de cada panel debería oscilar cerca de 1.80.

fig, axes = plt.subplots(1, 5, figsize=(14, 3.4), sharex=True, sharey=True)
xs = np.linspace(df1['metros'].min(), df1['metros'].max(), 60)

for j, (tr, te) in enumerate(k5.split(X1)):
    ax = axes[j]
    m = LinearRegression().fit(X1[tr], y1[tr])
    ax.scatter(X1[tr], y1[tr], c=COLOR_TRAIN, s=14, alpha=0.55)
    ax.scatter(X1[te], y1[te], c=COLOR_TEST, s=28, zorder=3)
    ax.plot(xs, m.predict(xs.reshape(-1, 1)), color=COLOR_LINEA, lw=1.3)
    ax.set_title(f'Fold {j+1}\npend = {m.coef_[0]:.2f}')
    ax.set_xlabel('m²')
    if j == 0:
        ax.set_ylabel('Precio')

fig.suptitle('Cinco rectas distintas: cada una ignora su bloque rojo', y=1.08)
plt.tight_layout()
plt.show()

Comparar RMSE y R² por fold

Barras del error y del \(R^2\) de cada bloque, y un scatter que cruza el RMSE con el número de atípicos en prueba: la hipótesis es que el fold que se los lleva tiene más error.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.1))

ax = axes[0]
ax.bar(tab1['fold'], tab1['RMSE'], color=COLOR_TRAIN, edgecolor='white')
ax.axhline(tab1['RMSE'].mean(), color=COLOR_TEST, ls='--',
           label=f'Promedio = {tab1["RMSE"].mean():.1f}')
ax.set_xlabel('Fold')
ax.set_ylabel('RMSE (miles de USD)')
ax.set_title('RMSE de prueba por fold')
ax.set_xticks(tab1['fold'])
ax.legend(frameon=False, fontsize=8)

ax = axes[1]
ax.bar(tab1['fold'], tab1['R2'], color='#54A24B', edgecolor='white')
ax.axhline(tab1['R2'].mean(), color=COLOR_TEST, ls='--',
           label=f'Promedio = {tab1["R2"].mean():.3f}')
ax.set_xlabel('Fold')
ax.set_ylabel('R²')
ax.set_title('R² de prueba por fold')
ax.set_xticks(tab1['fold'])
ax.legend(frameon=False, fontsize=8)

ax = axes[2]
ax.scatter(tab1['n_atipicos_en_prueba'], tab1['RMSE'],
           c=COLOR_TEST, s=70, zorder=3)
for _, r in tab1.iterrows():
    ax.annotate(f'F{int(r["fold"])}',
                (r['n_atipicos_en_prueba'] + 0.05, r['RMSE']),
                fontsize=8)
ax.set_xlabel('Atípicos en el fold de prueba')
ax.set_ylabel('RMSE')
ax.set_title('¿Tiene más error el fold que hereda atípicos?')
ax.set_xticks([0, 1, 2, 3])

plt.tight_layout()
plt.show()

Cómo leer la tabla y estas barras

No empieces por el promedio: empieza por el fold peor y el mejor, y busca una causa en las columnas de contexto (n_atipicos_en_prueba, rango de metros, desviación del precio en prueba).

  • El fold 1 se llevó los tres atípicos y su RMSE sube a 22.7, contra 14–16 en los otros cuatro. k-fold no está “siendo inestable porque sí”: te está avisando de sensibilidad a puntos influyentes.
  • La pendiente, aun así, se queda entre 1.76 y 1.86, pegada a la verdadera (1.80). La recta es estable; lo que se mueve es el error de predicción cuando los diamantes quedan fuera.
  • El \(R^2\) del fold 1 baja a 0.88 y el resto se queda en 0.91–0.96. Parte de esa caída es el RMSE más alto, no una población distinta de metros (el fold 1 cubre casi todo el rango, 46–164 m²).

El número que reportas es el promedio: RMSE 16.3 ± 3.2. El diagnóstico fold a fold es lo que te dice si te fías de ese promedio — y aquí el ±3.2 viene, sobre todo, de un solo bloque.

Ver predicciones y residuales cuando cada punto era prueba

Tres diagramas de dispersión out-of-fold (cada predicción se hizo cuando el punto estaba en el fold de prueba), coloreados por fold: predicho vs observado, residual vs predicho y residual vs metros. Cada punto fue evaluado por un modelo que no lo había visto.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.1))

ax = axes[0]
for j in range(1, 6):
    mask = fold_id == j
    ax.scatter(pred_cv[mask], y1[mask], c=colores[j - 1], s=28,
               label=f'Fold {j}', edgecolors='white', linewidths=0.4)
lo, hi = pred_cv.min(), pred_cv.max()
ax.plot([lo, hi], [lo, hi], color=COLOR_LINEA, lw=1.2)
ax.set_xlabel('Precio predicho (cuando el punto era prueba)')
ax.set_ylabel('Precio observado')
ax.set_title('Predicho vs observado out-of-fold')
ax.legend(frameon=False, fontsize=8)

ax = axes[1]
for j in range(1, 6):
    mask = fold_id == j
    ax.scatter(pred_cv[mask], resid_cv[mask], c=colores[j - 1], s=28,
               edgecolors='white', linewidths=0.4)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('Precio predicho out-of-fold')
ax.set_ylabel('Residual out-of-fold')
ax.set_title('Residuales de prueba, por fold')

ax = axes[2]
for j in range(1, 6):
    mask = fold_id == j
    ax.scatter(X1[mask, 0], resid_cv[mask], c=colores[j - 1], s=28,
               edgecolors='white', linewidths=0.4)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('Metros cuadrados')
ax.set_ylabel('Residual out-of-fold')
ax.set_title('¿El error crece con el tamaño?')

plt.tight_layout()
plt.show()

print('Residual out-of-fold de los atípicos:')
for i in idx_atipicos:
    print(f'  obs {i}: metros={X1[i,0]:.1f}, precio={y1[i]:.1f}, '
          f'pred_cv={pred_cv[i]:.1f}, residual={resid_cv[i]:.1f} '
          f'(fold {fold_id[i]})')

Residual out-of-fold de los atípicos:
  obs 18: metros=150.1, precio=336.4, pred_cv=301.6, residual=34.8 (fold 1)
  obs 61: metros=146.0, precio=237.2, pred_cv=294.1, residual=-56.8 (fold 1)
  obs 104: metros=47.9, precio=166.1, pred_cv=111.9, residual=54.2 (fold 1)

Estos tres scatters son la radiografía que un split único te esconde. Cada punto está pintado con el fold en el que fue prueba: su predicción la hizo un modelo que no lo había visto.

  • Los colores se mezclan alrededor de la diagonal: ningún bloque quedó sistemáticamente sesgado. Los tres puntos que más se alejan son los atípicos, todos del fold 1, con residuales de 34.8, −56.8 y 54.2: entre dos y tres veces y media el RMSE típico (~16). Eso explica el RMSE de 22.7 de ese fold.
  • Los residuales vs predicho no deberían dibujar un patrón en U (eso sería no linealidad). Un abanico a la derecha es la heterocedasticidad que ya vimos.
  • Residual vs metros: si la nube se abre en los m² altos, el RMSE global está siendo empujado por los apartamentos grandes. En un proyecto real, eso invita a ponderar, transformar o aceptar que el error no es constante.

Contrastar in-sample con k-fold (k=5 y k=10)

Se calcula el RMSE y el \(R^2\) de entrenamiento y se comparan con cross_val_score a 5 y a 10 folds. El in-sample es la referencia optimista; los dos CV deberían coincidir en el promedio.

rmse_in = np.sqrt(np.mean(resid_in**2))
rmse_cv5 = tab1['RMSE'].mean()

k10 = KFold(n_splits=10, shuffle=True, random_state=SEMILLA)
rmse_cv10 = -cross_val_score(
    LinearRegression(), X1, y1, cv=k10, scoring='neg_root_mean_squared_error'
)
r2_cv10 = cross_val_score(LinearRegression(), X1, y1, cv=k10, scoring='r2')
r2_cv5 = cross_val_score(LinearRegression(), X1, y1, cv=k5, scoring='r2')

print(f'RMSE in-sample (optimista):     {rmse_in:.2f}')
print(f'RMSE CV k=5:                    {rmse_cv5:.2f}  ± {tab1["RMSE"].std(ddof=0):.2f}')
print(f'RMSE CV k=10:                   {rmse_cv10.mean():.2f}  ± {rmse_cv10.std():.2f}')

# RMSE con los 120 residuales out-of-fold juntos (no promedio de RMSE por fold)
pred_cv10 = np.empty(n1)
for tr, te in k10.split(X1):
    pred_cv10[te] = LinearRegression().fit(X1[tr], y1[tr]).predict(X1[te])
print(f'RMSE CV k=10 (residuales juntos): {np.sqrt(np.mean((y1 - pred_cv10)**2)):.2f}')
print()
print(f'R2   in-sample:                 {ols1.rsquared:.3f}')
print(f'R2   CV k=5:                    {r2_cv5.mean():.3f}  ± {r2_cv5.std():.3f}')
print(f'R2   CV k=10:                   {r2_cv10.mean():.3f}  ± {r2_cv10.std():.3f}')
print()
print('R2 por fold (k=10):', np.round(r2_cv10, 3))
RMSE in-sample (optimista):     16.17
RMSE CV k=5:                    16.33  ± 3.24
RMSE CV k=10:                   16.08  ± 3.55
RMSE CV k=10 (residuales juntos): 16.47

R2   in-sample:                 0.939
R2   CV k=5:                    0.927  ± 0.031
R2   CV k=10:                   0.919  ± 0.045

R2 por fold (k=10): [0.84  0.897 0.931 0.883 0.957 0.961 0.951 0.957 0.963 0.852]

Comparar la dispersión del R² entre k=5 y k=10

Un boxplot por esquema: el promedio debería parecerse; lo que cambia es cuánto baila cada fold (con \(k=10\) cada prueba es más pequeña y más ruidosa).

fig, ax = plt.subplots(figsize=(7.5, 4.2))
ax.boxplot(
    [r2_cv5, r2_cv10],
    tick_labels=['k = 5', 'k = 10'],
    patch_artist=True,
    boxprops=dict(facecolor=COLOR_TRAIN, alpha=0.45),
    medianprops=dict(color=COLOR_TEST, lw=1.6),
)
ax.set_ylabel('R² de cada fold')
ax.set_title('Más folds: el promedio se mueve poco, el abanico se abre')
plt.tight_layout()
plt.show()

El RMSE in-sample (16.2) y el de CV con \(k=5\) (16.3) están casi pegados: el modelo es una recta, no tiene capacidad de memorizar 120 puntos, así que casi no hay sobreajuste. Lo que sí aparece es la dispersión entre folds (±3.2), empujada por el bloque de los atípicos. Con \(k=10\), el promedio de los RMSE por fold (16.08) sale incluso un poco por debajo del in-sample (16.17). No es que la validación cruzada se vuelva optimista: promediar los RMSE de los folds no es lo mismo que calcular un RMSE con todos los residuales de prueba juntos. Como la raíz cuadrada es cóncava, el promedio de las raíces queda por debajo de la raíz del promedio, y el efecto crece cuando los folds son pequeños (~12 puntos) y dispares. Con los 120 residuales out-of-fold juntos, el RMSE de \(k=10\) es 16.47, por encima del in-sample, como se espera. Los promedios de \(k=5\) y \(k=10\) casi coinciden (\(R^2\) de 0.927 y 0.919); lo que cambia es el abanico entre folds, que con \(k=10\) es más ancho (desviación de 0.045 contra 0.031, con un fold en 0.84), porque cada prueba es más pequeña y más ruidosa.

Pasa ahora a varias variables. Ahí k-fold no solo mueve el error: mueve todo el vector de coeficientes.

Ejemplo 2. Regresión multivariable: progresión de diabetes

El dataset diabetes de scikit-learn tiene 442 pacientes y 10 variables clínicas (edad, sexo, IMC, presión arterial y seis marcadores de sangre) para predecir un índice cuantitativo de progresión de la enfermedad un año después. Las columnas ya vienen centradas y escaladas —cada una tiene media 0 y suma de cuadrados 1; por eso su desviación estándar es 0.048 y no 1—: lo que ves no son kg/m² ni mmHg, sino desviaciones respecto al promedio en una escala artificial. Por eso los coeficientes salen en cientos: una unidad de esa escala equivale a unas 21 desviaciones estándar de la variable. El objetivo, en cambio, está en su escala original.

Con una sola variable podías mirar la recta. Con diez, k-fold se vuelve el instrumento para responder dos preguntas distintas: ¿cuánto error de predicción espero? y ¿los coeficientes se mantienen cuando cambio el bloque de pacientes que dejo fuera?

Cargar diabetes y mirar el resumen de las 10 variables

442 pacientes, 10 predictoras ya estandarizadas y un índice de progresión. Se imprime el describe y la correlación de cada variable con el objetivo.

datos = load_diabetes(as_frame=True)
X2 = datos.data.copy()
y2 = datos.target.copy()
df2 = X2.copy()
df2['progreso'] = y2

print(f'Observaciones: {len(df2)} | Predictoras: {X2.shape[1]}')
print('Variables:', list(X2.columns))
print()
print(df2.describe().round(3).T.to_string())
print()
print('Correlación de cada predictora con la progresión:')
print(X2.corrwith(y2).sort_values(ascending=False).round(3).to_string())
Observaciones: 442 | Predictoras: 10
Variables: ['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

          count     mean     std     min     25%      50%      75%      max
age       442.0   -0.000   0.048  -0.107  -0.037    0.005    0.038    0.111
sex       442.0    0.000   0.048  -0.045  -0.045   -0.045    0.051    0.051
bmi       442.0   -0.000   0.048  -0.090  -0.034   -0.007    0.031    0.171
bp        442.0   -0.000   0.048  -0.112  -0.037   -0.006    0.036    0.132
s1        442.0   -0.000   0.048  -0.127  -0.034   -0.004    0.028    0.154
s2        442.0    0.000   0.048  -0.116  -0.030   -0.004    0.030    0.199
s3        442.0   -0.000   0.048  -0.102  -0.035   -0.007    0.029    0.181
s4        442.0   -0.000   0.048  -0.076  -0.039   -0.003    0.034    0.185
s5        442.0    0.000   0.048  -0.126  -0.033   -0.002    0.032    0.134
s6        442.0    0.000   0.048  -0.138  -0.033   -0.001    0.028    0.136
progreso  442.0  152.133  77.093  25.000  87.000  140.500  211.500  346.000

Correlación de cada predictora con la progresión:
bmi    0.586
s5     0.566
bp     0.441
s4     0.430
s6     0.382
s1     0.212
age    0.188
s2     0.174
sex    0.043
s3    -0.395

Ver el objetivo, las correlaciones y dos predictoras fuertes

Histograma de la progresión, mapa de correlaciones (ahí se ve el bloque s1–s2) y un scatter bmi vs s5 coloreado por el progreso.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.0))

ax = axes[0]
ax.hist(y2, bins=20, color='#54A24B', edgecolor='white')
ax.set_xlabel('Progresión (índice clínico)')
ax.set_ylabel('Frecuencia')
ax.set_title('Distribución del objetivo')

ax = axes[1]
corr = df2.corr().values
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(len(df2.columns)))
ax.set_yticks(range(len(df2.columns)))
ax.set_xticklabels(df2.columns, rotation=90, fontsize=8)
ax.set_yticklabels(df2.columns, fontsize=8)
ax.set_title('Correlaciones')
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

ax = axes[2]
ax.scatter(X2['bmi'], X2['s5'], c=y2, s=22, cmap='viridis',
           edgecolors='white', linewidths=0.3)
ax.set_xlabel('bmi (estandarizado)')
ax.set_ylabel('s5 (estandarizado)')
ax.set_title('Dos predictoras fuertes, color = progreso')
cbar = fig.colorbar(ax.collections[0], ax=ax, fraction=0.046, pad=0.04)
cbar.set_label('progreso', fontsize=8)

plt.tight_layout()
plt.show()

Dibujar la relación marginal de bmi, s5 y bp con la progresión

Tres nubes, cada una con su recta simple. Ninguna variable sola va a explicar el fenómeno; el modelo multivariable suma esas señales débiles.

top = ['bmi', 's5', 'bp']
fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.0), sharey=True)
for ax, col in zip(axes, top):
    ax.scatter(X2[col], y2, c=COLOR_TRAIN, s=18, alpha=0.75,
               edgecolors='white', linewidths=0.3)
    m = LinearRegression().fit(X2[[col]].values, y2.values)
    xs = np.linspace(X2[col].min(), X2[col].max(), 50)
    ax.plot(xs, m.predict(xs.reshape(-1, 1)), color=COLOR_TEST, lw=1.5)
    r = X2[col].corr(y2)
    ax.set_xlabel(col)
    ax.set_title(f'{col}  (r = {r:.2f})')
axes[0].set_ylabel('Progresión')
fig.suptitle('Relación marginal de las tres predictoras más asociadas', y=1.03)
plt.tight_layout()
plt.show()

Cómo leer este dataset antes de modelar

El progreso cubre un rango amplio (de 25 a 346), es asimétrico a la derecha (media 152.1, mediana 140.5) y no está partido en dos nubes: tiene sentido tratarlo como regresión, no como clasificación disfrazada. En el mapa de correlaciones hay un bloque entre s1 y s2 (marcadores de sangre muy alineados): multicolinealidad. Eso no rompe las predicciones de una lineal, pero sí infla la varianza de esos coeficientes y, por tanto, los va a hacer bailar más de fold a fold. El scatter bmi vs s5 coloreado por el objetivo sugiere que las dos empujan en la misma dirección, sin un grupo raro que un fold pudiera acaparar.

Las nubes marginales de bmi, s5 y bp son lineales y ruidosas: ninguna variable sola va a entregar un \(R^2\) alto. El modelo multivariable existe precisamente para sumar esas señales débiles. k-fold dirá si esa suma se sostiene en pacientes que el modelo no vio.

Ajustar OLS multivariable con intercepto, muestra completa

Igual que en el ejemplo 1: add_constant pone la columna de unos y OLS estima los 10 coeficientes más el intercepto. El summary() es otra vez una referencia in-sample, no el veredicto.

X2_sm = sm.add_constant(X2)
ols2 = sm.OLS(y2, X2_sm).fit()
print(ols2.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 target   R-squared:                       0.518
Model:                            OLS   Adj. R-squared:                  0.507
Method:                 Least Squares   F-statistic:                     46.27
Date:                Sun, 27 Sep 2026   Prob (F-statistic):           3.83e-62
Time:                        20:55:52   Log-Likelihood:                -2386.0
No. Observations:                 442   AIC:                             4794.
Df Residuals:                     431   BIC:                             4839.
Df Model:                          10                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        152.1335      2.576     59.061      0.000     147.071     157.196
age          -10.0099     59.749     -0.168      0.867    -127.446     107.426
sex         -239.8156     61.222     -3.917      0.000    -360.147    -119.484
bmi          519.8459     66.533      7.813      0.000     389.076     650.616
bp           324.3846     65.422      4.958      0.000     195.799     452.970
s1          -792.1756    416.680     -1.901      0.058   -1611.153      26.802
s2           476.7390    339.030      1.406      0.160    -189.620    1143.098
s3           101.0433    212.531      0.475      0.635    -316.684     518.770
s4           177.0632    161.476      1.097      0.273    -140.315     494.441
s5           751.2737    171.900      4.370      0.000     413.407    1089.140
s6            67.6267     65.984      1.025      0.306     -62.064     197.318
==============================================================================
Omnibus:                        1.506   Durbin-Watson:                   2.029
Prob(Omnibus):                  0.471   Jarque-Bera (JB):                1.404
Skew:                           0.017   Prob(JB):                        0.496
Kurtosis:                       2.726   Cond. No.                         227.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Mira tres cosas del summary, no todas sus filas por igual. Primero, el \(R^2\) in-sample (0.518), que es un techo optimista. Segundo, qué coeficientes tienen un intervalo al 95% que no cruza cero: sex, bmi, bp y s5. Tercero, la multicolinealidad: aquí statsmodels no imprime su aviso (el número de condición es 227), pero el síntoma está en la tabla. s1 y s2, con correlación de 0.90 entre sí, salen con coeficientes grandes y de signo opuesto (−792 y 477) y, aun así, ninguno es significativo al 5% (p = 0.058 y 0.160): su error estándar está inflado (417 y 339, contra ~65 de bmi). k-fold va a mostrar cuánto se mueven esos coeficientes cuando cambias el bloque de pacientes que dejas fuera.

Diagnosticar el OLS in-sample y medirlo contra predecir la media

Predicho vs observado, residuales vs predicho y residuales vs IMC. Al final, el RMSE de “predecir siempre el promedio”: la línea base que k-fold tiene que ganar.

pred_in2 = ols2.fittedvalues
resid_in2 = ols2.resid

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.1))

ax = axes[0]
ax.scatter(pred_in2, y2, c=COLOR_TRAIN, s=18, alpha=0.75,
           edgecolors='white', linewidths=0.3)
lo, hi = pred_in2.min(), pred_in2.max()
ax.plot([lo, hi], [lo, hi], color=COLOR_LINEA, lw=1.2)
ax.set_xlabel('Progreso predicho')
ax.set_ylabel('Progreso observado')
ax.set_title('In-sample: predicho vs observado')

ax = axes[1]
ax.scatter(pred_in2, resid_in2, c=COLOR_TRAIN, s=18, alpha=0.75,
           edgecolors='white', linewidths=0.3)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('Progreso predicho')
ax.set_ylabel('Residual')
ax.set_title('In-sample: residuales')

ax = axes[2]
ax.scatter(X2['bmi'], resid_in2, c=COLOR_TRAIN, s=18, alpha=0.75,
           edgecolors='white', linewidths=0.3)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('bmi')
ax.set_ylabel('Residual')
ax.set_title('¿Queda patrón contra el IMC?')

plt.tight_layout()
plt.show()

print(f'R2 in-sample:   {ols2.rsquared:.3f}')
print(f'RMSE in-sample: {np.sqrt(np.mean(resid_in2**2)):.2f}')
print(f'RMSE de predecir la media: {np.sqrt(np.mean((y2 - y2.mean())**2)):.2f}')

R2 in-sample:   0.518
RMSE in-sample: 53.48
RMSE de predecir la media: 77.01

No hay una curva evidente en los residuales contra el IMC: la linealidad, en esa dirección, es razonable. El RMSE de “predecir siempre la media” es la línea base. Si en k-fold el modelo no le gana a esa línea de forma clara, el modelo multivariable no está pagando su complejidad.

Correr k-fold a mano: métricas, coeficientes y línea base por fold

Cinco rondas. En cada una se guarda MAE, RMSE, \(R^2\), el RMSE de predecir la media del entrenamiento, y el vector de 10 coeficientes. Esos coeficientes son los que después vamos a ver bailar o no.

k5b = KFold(n_splits=5, shuffle=True, random_state=SEMILLA)
fold_id2 = np.empty(len(y2), dtype=int)
pred_cv2 = np.empty(len(y2))
resid_cv2 = np.empty(len(y2))
coefs = []
filas2 = []

baseline_rmse = []

for j, (tr, te) in enumerate(k5b.split(X2)):
    fold_id2[te] = j + 1
    m = LinearRegression().fit(X2.iloc[tr], y2.iloc[tr])
    pred = m.predict(X2.iloc[te])
    pred_cv2[te] = pred
    resid_cv2[te] = y2.iloc[te].values - pred
    coefs.append(m.coef_)
    media_tr = y2.iloc[tr].mean()
    rmse_base = np.sqrt(mean_squared_error(y2.iloc[te], np.full(len(te), media_tr)))
    baseline_rmse.append(rmse_base)
    filas2.append({
        'fold': j + 1,
        'n_test': len(te),
        'progreso_media': y2.iloc[te].mean(),
        'progreso_sd': y2.iloc[te].std(ddof=0),
        'MAE': mean_absolute_error(y2.iloc[te], pred),
        'RMSE': np.sqrt(mean_squared_error(y2.iloc[te], pred)),
        'R2': r2_score(y2.iloc[te], pred),
        'RMSE_baseline': rmse_base,
    })

tab2 = pd.DataFrame(filas2)
coefs = pd.DataFrame(coefs, columns=X2.columns)
coefs.index = [f'fold {j}' for j in range(1, 6)]

print(tab2.round(2).to_string(index=False))
print()
print(f'RMSE CV: {tab2["RMSE"].mean():.2f}  ± {tab2["RMSE"].std(ddof=0):.2f}')
print(f'R2   CV: {tab2["R2"].mean():.3f}  ± {tab2["R2"].std(ddof=0):.3f}')
print(f'RMSE baseline (media del train): {tab2["RMSE_baseline"].mean():.2f}')
print(f'R2 in-sample: {ols2.rsquared:.3f}   |   R2 CV: {tab2["R2"].mean():.3f}')
 fold  n_test  progreso_media  progreso_sd   MAE  RMSE   R2  RMSE_baseline
    1      89          149.78        77.05 41.52 52.34 0.54          77.10
    2      89          143.76        73.45 43.59 53.47 0.47          74.20
    3      88          159.91        77.71 47.99 58.09 0.44          78.31
    4      88          148.93        78.94 43.24 53.91 0.53          79.05
    5      88          158.41        76.59 43.49 54.53 0.49          76.99

RMSE CV: 54.47  ± 1.95
R2   CV: 0.495  ± 0.037
RMSE baseline (media del train): 77.13
R2 in-sample: 0.518   |   R2 CV: 0.495

Comparar el modelo contra la línea base, fold a fold

Barras de RMSE (modelo vs media), barras de \(R^2\), y un scatter de \(R^2\) contra la desviación del progreso en ese fold: para no confundir un bloque “difícil” con uno que simplemente varía menos.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.1))

ax = axes[0]
ancho = 0.36
ax.bar(tab2['fold'] - ancho / 2, tab2['RMSE'], width=ancho,
       color=COLOR_TRAIN, label='Modelo lineal')
ax.bar(tab2['fold'] + ancho / 2, tab2['RMSE_baseline'], width=ancho,
       color='#B0B0B0', label='Baseline (media)')
ax.set_xlabel('Fold')
ax.set_ylabel('RMSE')
ax.set_title('El modelo vs predecir la media, fold a fold')
ax.set_xticks(tab2['fold'])
ax.legend(frameon=False, fontsize=8)

ax = axes[1]
ax.bar(tab2['fold'], tab2['R2'], color='#54A24B', edgecolor='white')
ax.axhline(tab2['R2'].mean(), color=COLOR_TEST, ls='--',
           label=f'Promedio = {tab2["R2"].mean():.3f}')
ax.set_xlabel('Fold')
ax.set_ylabel('R²')
ax.set_title('R² de prueba por fold')
ax.set_xticks(tab2['fold'])
ax.legend(frameon=False, fontsize=8)

ax = axes[2]
ax.scatter(tab2['progreso_sd'], tab2['R2'], c=COLOR_TEST, s=70, zorder=3)
for _, r in tab2.iterrows():
    ax.annotate(f'F{int(r["fold"])}',
                (r['progreso_sd'] + 0.4, r['R2']), fontsize=8)
ax.set_xlabel('Desviación del progreso en el fold de prueba')
ax.set_ylabel('R² de ese fold')
ax.set_title('Un fold “difícil” puede ser solo menos varianza en y')

plt.tight_layout()
plt.show()

Cómo leer el desempeño multivariable

El modelo le gana a la línea base en todos los folds (RMSE ~52–58 contra ~74–79). No es un empate afortunado de un solo split. Esa es la primera luz verde.

El \(R^2\) de CV (0.495) queda por debajo del in-sample (0.518). La brecha es el precio de generalizar; si fuera enorme, habría sobreajuste. Una lineal con 10 variables y 442 filas casi no tiene capacidad de memorizar, y la brecha es de apenas 0.02: señal de poco sobreajuste. Si el modelo es bueno o mediocre no lo dice esa brecha: lo dice la comparación con la línea base.

El scatter de la derecha es una trampa clásica de lectura: un fold con poca varianza en \(y\) tiene un \(SS_{tot}\) pequeño, y el \(R^2\) se vuelve inestable o bajo aunque el RMSE sea parecido. El fold 3, el de peor \(R^2\) (0.44), no tiene la desviación más baja; su RMSE sí es el más alto (58.1). Ahí el modelo predice peor de verdad, no es un artefacto del denominador.

Ver la tabla de coeficientes por fold

Diez columnas, cinco filas. Una columna estable es un efecto que sobrevive al cambio de pacientes; una que cambia de magnitud o de signo es un efecto frágil (o colineal).

print('Coeficientes por fold:')
print(coefs.round(1).to_string())
print()
print('Media ± sd entre folds:')
estabilidad = pd.DataFrame({
    'media': coefs.mean(),
    'sd': coefs.std(ddof=0),
    'min': coefs.min(),
    'max': coefs.max(),
    'ols_completo': ols2.params.drop('const').values,
})
print(estabilidad.round(1).to_string())
Coeficientes por fold:
         age    sex    bmi     bp     s1     s2     s3     s4     s5    s6
fold 1  -5.6 -253.9  507.6  347.7 -885.1  503.5   97.9  182.7  774.4  58.1
fold 2  17.6 -228.4  517.4  341.4 -693.1  329.8  113.7  242.6  709.6  85.5
fold 3 -31.9 -258.7  520.8  307.9 -694.5  432.0   81.9  199.4  733.6  38.8
fold 4 -25.4 -242.1  472.2  367.8 -915.7  624.7  179.4  167.9  788.4  94.9
fold 5  -6.7 -221.5  576.4  260.6 -732.1  444.8   27.9  112.9  730.6  64.7

Media ± sd entre folds:
     media    sd    min    max  ols_completo
age  -10.4  17.4  -31.9   17.6         -10.0
sex -240.9  14.3 -258.7 -221.5        -239.8
bmi  518.9  33.5  472.2  576.4         519.8
bp   325.1  37.6  260.6  367.8         324.4
s1  -784.1  96.5 -915.7 -693.1        -792.2
s2   466.9  96.7  329.8  624.7         476.7
s3   100.2  49.0   27.9  179.4         101.0
s4   181.1  42.3  112.9  242.6         177.1
s5   747.3  29.4  709.6  788.4         751.3
s6    68.4  20.0   38.8   94.9          67.6

Mapear esa estabilidad: heatmap y barras de error

A la izquierda, el mapa de color. A la derecha, la media entre folds ± su desviación, con el OLS de la muestra completa encima. Si la cruz cae dentro de la banda, el ajuste con toda la muestra es coherente con los de cada fold.

fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.6))

ax = axes[0]
im = ax.imshow(coefs.values, cmap='RdBu_r', aspect='auto')
ax.set_xticks(range(len(coefs.columns)))
ax.set_xticklabels(coefs.columns, rotation=90)
ax.set_yticks(range(5))
ax.set_yticklabels(coefs.index)
ax.set_title('Mapa de coeficientes (filas = folds)')
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

ax = axes[1]
xpos = np.arange(len(estabilidad))
ax.errorbar(
    xpos, estabilidad['media'], yerr=estabilidad['sd'],
    fmt='o', color=COLOR_TRAIN, ecolor=COLOR_TEST, capsize=3, ms=6,
    label='Media CV ± sd',
)
ax.scatter(xpos, estabilidad['ols_completo'], c=COLOR_LINEA, marker='x', s=40,
           label='OLS muestra completa', zorder=3)
ax.axhline(0, color='#888888', lw=0.8)
ax.set_xticks(xpos)
ax.set_xticklabels(estabilidad.index, rotation=90)
ax.set_ylabel('Coeficiente')
ax.set_title('Estabilidad: quién se mantiene y quién baila')
ax.legend(frameon=False, fontsize=8)

plt.tight_layout()
plt.show()

En el mapa, una columna de color uniforme es un efecto que sobrevive al cambio de pacientes: bmi (media 519, sd 34), s5 (747, sd 29) y sex (−241, sd 14). Una columna que cambia de intensidad —s1 y s2, con sd cercana a 97— es el síntoma de dos marcadores que se cuentan la misma historia: el OLS de la muestra completa les pone coeficientes grandes y de signo opuesto. k-fold muestra que el signo opuesto sí se repite en los cinco bloques (s1 entre −916 y −693; s2 entre 330 y 625): lo que se mueve es la magnitud. La única variable que cambia de signo entre folds es age (de −31.9 a 17.6), coherente con su intervalo del summary, que cruza cero.

El panel de la derecha resume lo mismo con barras de error. Las cruces del OLS completo caen dentro de la banda de k-fold, pero eso era de esperar: los folds comparten el 75% de sus filas de entrenamiento, así que la desviación entre folds subestima la incertidumbre real de cada coeficiente (para s1, 96.5 entre folds contra un error estándar de 417 en el summary). La estabilidad entre folds es una condición necesaria para interpretar un coeficiente, no una prueba de que el efecto exista. La decisión práctica se parte en dos. Para predecir, el modelo sirve: RMSE de CV ~54 contra baseline ~77. Para interpretar, quédate con bmi, bp, s5 y sex; no leas un relato clínico sobre s1 y s2 como si fueran efectos separados.

Ver predicciones out-of-fold y si algún fold se llevó otra población

Predicho vs observado (ahora sí, honesto), residuales de prueba, y boxplots del progreso observado en cada bloque. Los boxplots deberían parecerse: si un fold se llevó solo a los pacientes graves, su RMSE no es comparable.

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.1))

ax = axes[0]
for j in range(1, 6):
    mask = fold_id2 == j
    ax.scatter(pred_cv2[mask], y2.values[mask], c=colores[j - 1], s=18,
               label=f'Fold {j}', edgecolors='white', linewidths=0.3, alpha=0.9)
lo, hi = pred_cv2.min(), pred_cv2.max()
ax.plot([lo, hi], [lo, hi], color=COLOR_LINEA, lw=1.2)
ax.set_xlabel('Predicho out-of-fold')
ax.set_ylabel('Observado')
ax.set_title('Predicho vs observado (k-fold)')
ax.legend(frameon=False, fontsize=8)

ax = axes[1]
for j in range(1, 6):
    mask = fold_id2 == j
    ax.scatter(pred_cv2[mask], resid_cv2[mask], c=colores[j - 1], s=18,
               edgecolors='white', linewidths=0.3, alpha=0.9)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('Predicho out-of-fold')
ax.set_ylabel('Residual out-of-fold')
ax.set_title('Residuales de prueba')

ax = axes[2]
partes = [y2.values[fold_id2 == j] for j in range(1, 6)]
bp = ax.boxplot(partes, tick_labels=[str(j) for j in range(1, 6)],
                patch_artist=True, medianprops=dict(color=COLOR_TEST, lw=1.4))
for caja in bp['boxes']:
    caja.set_facecolor(COLOR_TRAIN)
    caja.set_alpha(0.35)
ax.set_xlabel('Fold de prueba')
ax.set_ylabel('Progreso observado')
ax.set_title('¿Algún fold se llevó una población distinta?')

plt.tight_layout()
plt.show()

Localizar dónde se equivoca el modelo

Nube bmi–progreso coloreada por el residual CV, residual vs s5, y los cinco pacientes con mayor |residual| cuando fueron prueba. El RMSE es un promedio; aquí se ven las colas.

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))

ax = axes[0]
ax.scatter(X2['bmi'], y2, c=resid_cv2, s=22, cmap='coolwarm',
           vmin=-np.quantile(np.abs(resid_cv2), 0.95),
           vmax=np.quantile(np.abs(resid_cv2), 0.95),
           edgecolors='white', linewidths=0.3)
ax.set_xlabel('bmi')
ax.set_ylabel('Progreso')
ax.set_title('Nube bmi–progreso, color = residual CV')
cbar = fig.colorbar(ax.collections[0], ax=ax, fraction=0.046, pad=0.04)
cbar.set_label('residual out-of-fold', fontsize=8)

ax = axes[1]
ax.scatter(X2['s5'], resid_cv2, c=COLOR_TRAIN, s=18, alpha=0.75,
           edgecolors='white', linewidths=0.3)
ax.axhline(0, color=COLOR_LINEA, lw=1.0)
ax.set_xlabel('s5')
ax.set_ylabel('Residual out-of-fold')
ax.set_title('¿Queda señal en s5 después del modelo?')

plt.tight_layout()
plt.show()

abs_res = pd.Series(np.abs(resid_cv2), index=X2.index, name='|residual CV|')
peores = abs_res.nlargest(5)
print('Cinco pacientes con mayor |residual| cuando fueron prueba:')
print(
    df2.loc[peores.index, ['bmi', 'bp', 's5', 'progreso']]
    .assign(residual_cv=resid_cv2[peores.index], fold=fold_id2[peores.index])
    .round(3)
    .to_string()
)

Cinco pacientes con mayor |residual| cuando fueron prueba:
       bmi     bp     s5  progreso  residual_cv  fold
56   0.041 -0.026  0.034      52.0     -158.315     1
102  0.028 -0.016 -0.005     302.0      153.959     2
92   0.005  0.098  0.013      48.0     -150.868     4
304 -0.039  0.053 -0.023     253.0      136.885     5
359  0.006  0.032  0.015     311.0      131.012     2

El predicho vs observado out-of-fold es más honesto que el in-sample: la nube se abre un poco más: es la brecha entre el \(R^2\) in-sample (0.518) y el de CV (0.495). Los boxplots del progreso por fold se parecen (medias entre 143.8 y 159.9): ningún bloque se llevó solo a los pacientes graves, así que sus RMSE son comparables.

Los dos scatters de la celda anterior cierran el diagnóstico:

  • Colorear la nube bmi–progreso por el residual CV te dice dónde se equivoca el modelo. El degradado vertical (rojo arriba, azul abajo) es esperable por construcción: el residual es \(y - \hat{y}\), así que los valores altos de \(y\) tienden a quedar subestimados. Lo que delataría un patrón no capturado es una franja de un solo color a lo largo de bmi (no linealidad, un subgrupo).
  • Residual vs s5: una pendiente lineal casi no puede quedar, porque el modelo ya incluye a s5; lo que habría que buscar es curvatura, que pediría un término no lineal, como los splines del próximo artículo de la serie. Aquí la nube está centrada en cero y sin curva visible: s5 ya hizo su trabajo.

Los cinco pacientes con mayor |residual| out-of-fold no son “errores de código”: son los casos que un modelo lineal promedio no explica. En un proyecto, los abrirías uno a uno (dato mal cargado, comorbilidad, un perfil clínico distinto). Aquí sirven para recordar que el RMSE de CV es un promedio que esconde colas: el peor residual (−158.3) es casi tres veces el RMSE de CV (54.5).

Confirmar con scikit-learn: in-sample, k=5, k=10 y la baseline

Una tabla con \(R^2\) y RMSE de los cuatro esquemas, más el vector de \(R^2\) de cada fold. Es el mismo k-fold de arriba, ahora con cross_val_score.

r2_cv5_b = cross_val_score(LinearRegression(), X2, y2, cv=k5b, scoring='r2')
k10b = KFold(n_splits=10, shuffle=True, random_state=SEMILLA)
r2_cv10_b = cross_val_score(LinearRegression(), X2, y2, cv=k10b, scoring='r2')
rmse_cv5_b = -cross_val_score(
    LinearRegression(), X2, y2, cv=k5b, scoring='neg_root_mean_squared_error'
)
rmse_cv10_b = -cross_val_score(
    LinearRegression(), X2, y2, cv=k10b, scoring='neg_root_mean_squared_error'
)

comp = pd.DataFrame({
    'esquema': ['In-sample', 'CV k=5', 'CV k=10', 'Baseline (media)'],
    'R2': [
        ols2.rsquared,
        r2_cv5_b.mean(),
        r2_cv10_b.mean(),
        np.nan,
    ],
    'RMSE': [
        np.sqrt(np.mean(resid_in2**2)),
        rmse_cv5_b.mean(),
        rmse_cv10_b.mean(),
        tab2['RMSE_baseline'].mean(),
    ],
})
print(comp.round(3).to_string(index=False))
print()
print('R2 k=5 :', np.round(r2_cv5_b, 3), f'→ {r2_cv5_b.mean():.3f} ± {r2_cv5_b.std():.3f}')
print('R2 k=10:', np.round(r2_cv10_b, 3), f'→ {r2_cv10_b.mean():.3f} ± {r2_cv10_b.std():.3f}')
         esquema    R2   RMSE
       In-sample 0.518 53.476
          CV k=5 0.495 54.467
         CV k=10 0.490 54.391
Baseline (media)   NaN 77.130

R2 k=5 : [0.539 0.47  0.441 0.534 0.493] → 0.495 ± 0.037
R2 k=10: [0.543 0.507 0.47  0.479 0.551 0.344 0.534 0.525 0.503 0.445] → 0.490 ± 0.058

Comparar el abanico del R² con k=5 y con k=10

Boxplot final del ejemplo multivariable. Con \(n=442\) el promedio debería coincidir; el abanico de \(k=10\) suele ser más ancho porque cada prueba es más pequeña.

fig, ax = plt.subplots(figsize=(7.5, 4.2))
ax.boxplot(
    [r2_cv5_b, r2_cv10_b],
    tick_labels=['k = 5', 'k = 10'],
    patch_artist=True,
    boxprops=dict(facecolor=COLOR_TRAIN, alpha=0.45),
    medianprops=dict(color=COLOR_TEST, lw=1.6),
)
ax.set_ylabel('R² de cada fold')
ax.set_title('Diabetes: k=5 y k=10 cuentan la misma historia')
plt.tight_layout()
plt.show()

In-sample gana en \(R^2\) (0.518) y pierde en honestidad. \(k=5\) y \(k=10\) coinciden en el promedio (~0.49). El abanico de \(k=10\) es más ancho (aparece un fold en 0.34), no más fino: cada prueba tiene ~44 pacientes y el \(R^2\) de un bloque pequeño es más ruidoso. Ese abanico describe folds individuales, no la precisión del promedio: los dos promedios son defendibles. Aquí reportarías el de \(k=5\) (0.495) y usarías el de \(k=10\) (0.490) para confirmar que el resultado no dependía de haber cortado en cinco.

Y ese 0.49 solo es interesante porque gana a predecir la media (RMSE 54 contra 77). En un fenómeno biológico ruidoso, un \(R^2\) de CV cercano a 0.5 es un modelo útil, no uno mediocre.

Cierre técnico

k-fold no es un trámite de scikit-learn. Es un diseño para medir el error de predicción sin apostar a una sola partición. Lo que te tiene que quedar:

  • Dibuja las particiones al menos una vez. Si no puedes señalar quién entrena y quién prueba en cada ronda, todavía no entendiste el método.
  • Baraja cuando las filas sean intercambiables. Si los datos vienen ordenados, k-fold sin shuffle mide otro experimento.
  • Lee el fold peor, no solo el promedio. Casi siempre hay una causa: un atípico, un rango distinto de \(x\), menos varianza en \(y\).
  • Reporta RMSE y \(R^2\) juntos, y compara contra una baseline (predecir la media). Sin esa ancla, 0.49 no dice nada.
  • En regresiones multivariables, mira también la estabilidad de los coeficientes. Predecir bien y poder interpretar no son el mismo premio.
  • \(k=5\) o \(k=10\) basta en \(n\) de unos cientos. No uses la misma validación cruzada para elegir el modelo y para reportar su desempeño.

En el artículo anterior, de evaluación de modelos y bootstrap, vimos el panorama completo (métricas, una partición, validación cruzada y bootstrap). Este se queda en una sola herramienta, a propósito, para que la puedas ver por dentro. En el próximo artículo, sobre splines, k-fold pasa a su uso más frecuente: elegir, solo con el entrenamiento, cuánta flexibilidad darle a un modelo.

💬 ¿Te sirvió?

Deja en los comentarios una duda o un caso donde aplicarías esto — respondo todos. Sígueme para no perderte el próximo artículo de la serie y comparte con alguien que esté aprendiendo análisis de datos.

👉 El código completo está disponible para ejecutar directamente.