4.5 Regresión polinomial

Las relaciones lineales no siempre son suficientes. La regresión polinomial añade potencias de un predictor al modelo, permitiendo que el resultado esperado se curve. Esta sección discute cómo elegir el grado del polinomio y por qué demasiada flexibilidad lleva al sobreajuste.

4.5.1 El modelo polinomial

El modelo de regresión polinomial añade términos cuadrados, cúbicos o de orden superior:

$\mu_i = \alpha + \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3$

El término cuadrático ( $\beta_2 x_i^2$ ) permite que la relación se curve. El término cúbico ( $\beta_3 x_i^3$ ) permite un punto de inflexión. Los términos de orden superior añaden más flexibilidad pero también más riesgo de sobreajuste.

# Polynomial regression: height vs weight
# The relationship is slightly nonlinear for the full age range

set.seed(42)
n <- 200
weight <- runif(n, 10, 65)

# True relationship: nonlinear (growth curve)

height <- 50 + 3 _ weight - 0.02 _ weight^2 + rnorm(n, sd = 3)

# Fit linear, quadratic, and cubic models

m1 <- lm(height ~ weight)
m2 <- lm(height ~ weight + I(weight^2))
m3 <- lm(height ~ weight + I(weight^2) + I(weight^3))

# Compare R-squared

cat("Linear: R-squared =", round(summary(m1)$r.squared, 4), "\n")
cat("Quadratic: R-squared =", round(summary(m2)$r.squared, 4), "\n")
cat("Cubic: R-squared =", round(summary(m3)$r.squared, 4), "\n")

# Plot

plot(weight, height, pch = 16, col = "steelblue",
main = "Polynomial regression: linear vs quadratic vs cubic",
xlab = "Weight (kg)", ylab = "Height (cm)")
x_seq <- seq(min(weight), max(weight), length.out = 100)
lines(x_seq, predict(m1, data.frame(weight = x_seq)), lwd = 2, col = "firebrick")
lines(x_seq, predict(m2, data.frame(weight = x_seq)), lwd = 2, col = "darkgreen")
lines(x_seq, predict(m3, data.frame(weight = x_seq)), lwd = 2, col = "purple")
legend("bottomright",
legend = c("Linear", "Quadratic", "Cubic"),
col = c("firebrick", "darkgreen", "purple"), lwd = 2)

Click Run to execute this code in your browser with WebR

4.5.2 Estandarización del predictor

Al usar términos polinomiales, los valores del predictor pueden volverse muy grandes (e.g., $65^3 = 274625$ ). Esto causa problemas numéricos. La solución es estandarizar el predictor antes de calcular las potencias:

$z_i = \frac{x_i - \bar{x}}{s_x}$

Entonces el modelo se convierte en:

$\mu_i = \alpha + \beta_1 z_i + \beta_2 z_i^2 + \beta_3 z_i^3$

# Standardize the predictor before polynomial regression
set.seed(42)
n <- 200
weight <- runif(n, 10, 65)
height <- 50 + 3 * weight - 0.02 * weight^2 + rnorm(n, sd = 3)

# Standardize weight

weight_z <- (weight - mean(weight)) / sd(weight)

# Fit models with standardized predictor

m1z <- lm(height ~ weight_z)
m2z <- lm(height ~ weight_z + I(weight_z^2))
m3z <- lm(height ~ weight_z + I(weight_z^2) + I(weight_z^3))

cat("With standardized predictor:\n")
cat("Linear: R-squared =", round(summary(m1z)$r.squared, 4), "\n")
cat("Quadratic: R-squared =", round(summary(m2z)$r.squared, 4), "\n")
cat("Cubic: R-squared =", round(summary(m3z)$r.squared, 4), "\n")

# Coefficients are now on a comparable scale

cat("\nQuadratic coefficients:\n")
cat("alpha:", round(coef(m2z)[1], 2), "\n")
cat("beta1:", round(coef(m2z)[2], 2), "\n")
cat("beta2:", round(coef(m2z)[3], 2), "\n")

Click Run to execute this code in your browser with WebR

4.5.3 El peligro del sobreajuste

Los polinomios de orden superior pueden oscilar salvajemente entre los puntos de datos, ajustando ruido en lugar de señal. Esto es sobreajuste — el modelo ajusta bien la muestra pero generaliza pobremente.

# Demonstrate overfitting with high-degree polynomials
set.seed(42)
n <- 30
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)

# Fit polynomials of increasing degree

par(mfrow = c(2, 2))
for (deg in c(1, 3, 5, 10)) {
fit <- lm(y ~ poly(x, degree = deg, raw = TRUE))

plot(x, y, pch = 16, col = "steelblue",
main = sprintf("Degree %d (R² = %.3f)", deg, summary(fit)$r.squared),
xlab = "x", ylab = "y", ylim = c(-5, 25))
x_seq <- seq(0, 10, length.out = 200)
lines(x_seq, predict(fit, data.frame(x = x_seq)), lwd = 2, col = "firebrick")
}
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

El polinomio de grado 10 ajusta los datos casi perfectamente, pero oscila salvajemente entre los puntos de datos. Este modelo generalizaría terriblemente a datos nuevos. El modelo de grado 1 (lineal) captura la relación verdadera sin sobreajustar.

Rethinking

La simplicidad es una virtud. El gólem no sabe cuándo dejar de añadir parámetros. Ajustará felizmente un polinomio de grado 20 si se lo pide. Es trabajo del modelador elegir el modelo más simple que capture las características importantes de los datos. Este es el principio de parsimonia: usar la menor complejidad necesaria. Los criterios de información (Capítulo 6) proporcionan una forma principiada de equilibrar ajuste y complejidad.

Overthinking

Splines y polinomios naturales. La regresión polinomial es la forma más simple de regresión no lineal, pero tiene limitaciones: los polinomios de alto grado son inestables en los bordes de los datos. Los splines son polinomios por tramos que son más suaves y estables. Un spline divide el rango del predictor en segmentos y ajusta un polinomio de bajo grado a cada segmento, uniéndolos suavemente en nudos. El paquete rethinking proporciona spline() y funciones relacionadas. Los splines son más flexibles que los polinomios y menos propensos a efectos de borde, pero requieren elegir el número y la ubicación de los nudos.

Resumen de la Sección 4.5

La regresión polinomial añade potencias del predictor: $\mu_i = \alpha + \beta_1 x_i + \beta_2 x_i^2 + \ldots$
Estandarice el predictor antes de calcular potencias para evitar problemas numéricos.
Los polinomios de alto grado sobreajustan: oscilan salvajemente entre los puntos de datos.
Use el modelo más simple que capture las características importantes — la parsimonia es una virtud.
Los splines son una alternativa más estable a los polinomios de alto grado.

Añadiendo un predictor

Resumen