library(ISLR2)
library(splines)
library(mgcv)
library(car)
library(strucchange)
library(ggplot2)
data(Wage)
data("mcycle", package = "MASS")
cv_mse <- function(formula, data, K = 10, seed = 1) {
set.seed(seed)
folds <- sample(rep(1:K, length.out = nrow(data)))
y <- eval(formula[[2]], data)
errs <- sapply(1:K, function(k) {
fit <- lm(formula, data = data[folds != k, ])
pred <- predict(fit, newdata = data[folds == k, ])
mean((y[folds == k] - pred)^2)
})
mean(errs)
}
pred_band <- function(fit, newdata) {
p <- predict(fit, newdata, se.fit = TRUE)
data.frame(newdata, fit = p$fit,
lwr = p$fit - 2 * p$se.fit, upr = p$fit + 2 * p$se.fit)
}
plot_band <- function(df, titulo, color = "darkred") {
ggplot(df, aes(age, fit)) +
geom_point(data = Wage, aes(age, wage), inherit.aes = FALSE,
alpha = 0.15, color = "grey40") +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.25, fill = color) +
geom_line(color = color, linewidth = 1) +
labs(title = titulo, x = "age", y = "wage") +
theme_minimal()
}
age_grid <- data.frame(age = seq(min(Wage$age), max(Wage$age), length.out = 200))Regresión no lineal
Polinomios, escalones, cambio estructural, splines, regresión local y GAMs
Ajuste con banda de confianza
plot_band(pred_band(fits[[4]], age_grid), "Polinomio de grado 4")
Cuidado: colas y estructura local
fit_m10 <- lm(accel ~ poly(times, 10), data = mcycle)
tt <- data.frame(times = seq(min(mcycle$times), max(mcycle$times), length.out = 300))
plot(mcycle$times, mcycle$accel, pch = 20, col = "grey40", xlab = "times (ms)", ylab = "accel")
lines(tt$times, predict(fit_m10, tt), col = "red", lwd = 2)
lines(tt$times, predict(loess(accel ~ times, mcycle, span = 0.3), tt), col = "blue", lwd = 2)
legend("bottomright", c("Polinomio gr. 10", "loess"), col = c("red", "blue"), lwd = 2)
Los polinomios son globales: un dato en un extremo mueve toda la curva y las colas oscilan.
Cortes elegidos por el analista
cortes <- c(18, 30, 45, 60, 80)
fit_step2 <- lm(wage ~ cut(age, breaks = cortes, include.lowest = TRUE), data = Wage)
head(model.matrix(fit_step2), 3) (Intercept) cut(age, breaks = cortes, include.lowest = TRUE)(30,45]
231655 1 0
86582 1 0
161300 1 1
cut(age, breaks = cortes, include.lowest = TRUE)(45,60]
231655 0
86582 0
161300 0
cut(age, breaks = cortes, include.lowest = TRUE)(60,80]
231655 0
86582 0
161300 0plot_band(pred_band(fit_step2, age_grid), "Función escalón (cortes en 30, 45, 60)", color = "darkgreen")
Cambio estructural con una variable binaria
\[ y_t = \beta_0 + \beta_1 t + \gamma_0 D_t + \gamma_1\, t D_t + \varepsilon_t, \qquad D_t = 1\{t > T_0\} \]
- \(\gamma_0\): salto en el intercepto en \(T_0\) (si centramos \(t\) en \(T_0\))
- \(\gamma_1\): cambio en la pendiente (la “deriva” de la tendencia)
Datos: caudal anual del Nilo (1871–1970), presa de Asuán en 1898.
nile <- data.frame(year = as.numeric(time(Nile)), flow = as.numeric(Nile))
T0 <- 1898
nile$t <- nile$year - T0
nile$D <- as.numeric(nile$year > T0)
plot(nile$year, nile$flow, type = "l", xlab = "año", ylab = "caudal")
abline(v = T0, lty = 2)
¿Y si no conocemos la fecha del quiebre?
Se calcula el estadístico \(F\) para cada fecha posible y se toma el supremo (sup-F). Su distribución no es F: usar valores críticos de strucchange.
fs <- Fstats(Nile ~ 1)
sctest(fs)
supF test
data: fs
sup.F = 75.93, p-value = 2.22e-16plot(fs)
breakdates(breakpoints(Nile ~ 1, breaks = 1))[1] 1898Bandas de confianza
plot_band(pred_band(fit_bs, age_grid), "B-spline (nudos 25, 40, 60)")
plot_band(pred_band(fit_ns, age_grid), "Spline natural (4 gl)", color = "steelblue")
Spline de suavizado en R
ss_gcv <- smooth.spline(Wage$age, Wage$wage) # lambda por GCV
c(df = ss_gcv$df, lambda = ss_gcv$lambda) df lambda
6.4682992 0.0348627 dfs <- c(2, 5, 10, 30)
ss_list <- lapply(dfs, function(g) smooth.spline(Wage$age, Wage$wage, df = g))
plot(Wage$age, Wage$wage, pch = 20, col = gray(0.6, 0.25), xlab = "age", ylab = "wage")
for (i in seq_along(dfs)) lines(ss_list[[i]], col = i + 1, lwd = 2)
legend("topright", paste("gl =", dfs), col = 2:5, lwd = 2)
Todos los métodos sobre Wage
modelos <- list("Polinomio grado 4" = fits[[4]],
"B-spline (3 nudos)" = fit_bs,
"Spline natural (4 gl)" = fit_ns)
curvas <- do.call(rbind, lapply(names(modelos), function(nm)
data.frame(modelo = nm, age = age_grid$age, fit = predict(modelos[[nm]], age_grid))))
curvas <- rbind(curvas,
data.frame(modelo = "Spline de suavizado (GCV)", age = age_grid$age,
fit = predict(ss_gcv, age_grid$age)$y))
ggplot(curvas, aes(age, fit, color = modelo)) +
geom_line(linewidth = 1) +
labs(x = "age", y = "wage", color = NULL) +
theme_minimal() + theme(legend.position = "bottom")
loess() en R
lo_default <- loess(accel ~ times, data = mcycle) # span = 0.75, grado 2
lo_03 <- loess(accel ~ times, data = mcycle, span = 0.3)
plot(mcycle$times, mcycle$accel, pch = 20, col = "grey40", xlab = "times", ylab = "accel")
lines(tt$times, predict(lo_default, tt), col = "red", lwd = 2)
lines(tt$times, predict(lo_03, tt), col = "blue", lwd = 2)
legend("bottomright", c("span = 0.75", "span = 0.3"), col = c("red", "blue"), lwd = 2)
Comparación con loess()
x <- mcycle$times; y <- mcycle$accel
grid <- seq(min(x), max(x), length.out = 200)
mi_loc <- sapply(grid, local_fit, x = x, y = y, span = 0.3, degree = 2)
lo_dir <- loess(accel ~ times, data = mcycle, span = 0.3, degree = 2,
control = loess.control(surface = "direct"))
lo_pred <- predict(lo_dir, data.frame(times = grid))
max(abs(mi_loc - lo_pred))[1] 1.847411e-13plot(x, y, pch = 20, col = "grey40", xlab = "times", ylab = "accel")
lines(grid, mi_loc, col = "red", lwd = 3)
lines(grid, lo_pred, col = "blue", lwd = 1, lty = 2)
legend("bottomright", c("Mi implementación", "loess"), col = c("red", "blue"), lwd = c(3, 1), lty = c(1, 2))
Debe ser prácticamente cero. Si queda una diferencia mínima, se debe a detalles de implementación en la elección del \(k\)-ésimo vecino. Con surface = "interpolate" (el valor por defecto de loess) se interpola sobre una malla y la diferencia es mayor.
Efecto del span y del grado
par(mfrow = c(1, 2))
for (gr in 1:2) {
plot(x, y, pch = 20, col = "grey60", main = paste("Grado", gr), xlab = "times", ylab = "accel")
for (i in 1:3) {
sp <- c(0.1, 0.3, 0.6)[i]
lines(grid, sapply(grid, local_fit, x = x, y = y, span = sp, degree = gr), col = i + 1, lwd = 2)
}
legend("bottomright", paste("span =", c(0.1, 0.3, 0.6)), col = 2:4, lwd = 2, cex = 0.8)
}par(mfrow = c(1, 1))
span chico: poco sesgo, mucha varianza. span grande: lo contrario. Es el mismo dilema que df en splines.
Elegir el span por LOOCV
loocv_span <- function(span, degree = 2) {
mean(sapply(seq_along(y), function(i)
(y[i] - local_fit(x[i], x[-i], y[-i], span, degree))^2))
}
spans <- seq(0.15, 0.8, by = 0.05)
cv_span <- sapply(spans, loocv_span)
plot(spans, cv_span, type = "b", xlab = "span", ylab = "LOOCV")spans[which.min(cv_span)][1] 0.35
Grados de libertad efectivos de loess: traza de la matriz suavizadora.
lo_03$trace.hat[1] 12.43856Los efectos estimados
plot(gam1, pages = 1, all.terms = TRUE, shade = TRUE)
¿Recupera el GAM la forma verdadera?
malla <- data.frame(edad = seq(18, 80, length.out = 200), expo = 1)
malla$estimada <- predict(gam_p, malla, type = "response")
malla$verdadera <- freq_verdadera(malla$edad)
ggplot(malla, aes(edad)) +
geom_line(aes(y = estimada, color = "GAM estimado"), linewidth = 1) +
geom_line(aes(y = verdadera, color = "Verdadera"), linewidth = 1, linetype = 2) +
labs(y = "Frecuencia esperada (exposición = 1)", color = NULL) +
theme_minimal()
Sin imponer una forma cuadrática, el GAM aprende la U. Es el análogo flexible del turning point que probamos con el método delta.