Métodos de regresión no lineal: regresión polinómica
Introducción
Los modelos lineales tienen la ventaja de ser fácilmente interpretables, sin embargo, pueden tener limitaciones importantes en capacidad predictiva. Esto se debe a que, la asunción de linealidad, es con frecuencia una aproximación demasiado simple para describir las relaciones reales entre variables. A continuación, se describen métodos que permiten relajar la condición de linealidad intentando mantener al mismo tiempo una interpretabilidad alta.
Regresión polinómica: Consigue añadir curvatura al modelo introduciendo nuevos predictores que se obtienen al elevar todos o algunos de los predictores originales a distintas potencias.
Regresión no lineal con un único predictor
Regresión polinómica
La forma más sencilla de incorporar flexibilidad a un modelo lineal es introduciendo nuevos predictores obtenidos al elevar a distintas potencias el predictor original.
Partiendo del modelo lineal
\[ y_i = \beta_0 + \beta_1x_i + \epsilon_i \]
Se obtiene un modelo polinómico de grado d a partir de la ecuación
\[ y_i = \beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i+ \epsilon_i \]
Los modelos polinómicos se pueden ajustar mediante regresión lineal por mínimos cuadrados ya que, aunque generan modelos no lineales, su ecuación no deja de ser una ecuación lineal con predictores x, \[x^2, x^3, ..., x^d.\] Por esta misma razón, las funciones polinómicas pueden emplearse en regresión logística para predecir respuestas binarias. Solo es necesario realizar una transformación logit.
\[ P(y_i>Y|x_i=X) = \frac{exp(\beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i)}{1 + exp(\beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i)} \]
En el libro Introduction to Statistical Learning desaconsejan el uso de modelos polinómicos con grado mayor de 3 o 4 debido a un exceso de flexibilidad (overfitting), principalmente en los extremos del predictor X. La selección del grado de polinomio óptimo puede hacerse mediante cross validation.
Ejemplo 1
Los datos contienen información de las entidades del país. Considerando ocho variables, entre las cuales se encuentra la producción de carne en toneladas.
setwd("~/EALMV9/U2/U2A5")
library(pacman)## Warning: package 'pacman' was built under R version 4.0.5library(readxl)
p_load("base64enc", "htmltools", "mime", "xfun", "prettydoc","readr", "knitr","DT","dplyr", "ggplot2")Descarga de datos
## Disponible para descargar el código
xfun::embed_file("U2A5.Rmd")## Datos disponibles datos utilizados en la actividad
xfun::embed_file("bovinoNacional2019.xlsx")Download bovinoNacional2019.xlsx
datos <- read_excel("bovinoNacional2019.xlsx")
datos$Precipitaciones[1:5]## [1] 477 199 184 1202 1989poly(datos$Precipitaciones, degree = 4, raw = TRUE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 477 227529 108531333 5.176945e+10
## [2,] 199 39601 7880599 1.568239e+09
## [3,] 184 33856 6229504 1.146229e+09
## [4,] 1202 1444804 1736654408 2.087459e+12
## [5,] 1989 3956121 7868724669 1.565089e+13poly(datos$Precipitaciones, degree = 4, raw = FALSE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] -0.1531461 0.08889459 0.05645864 -0.1574215
## [2,] -0.2565531 0.32787277 -0.34413098 0.2638193
## [3,] -0.2621326 0.34305199 -0.37716497 0.3152933
## [4,] 0.1165305 -0.15632161 -0.10639730 0.1497341
## [5,] 0.4092690 0.19607835 -0.29731475 -0.6223532modelo_poli4 <- lm(ProduccionCarne ~ poly(Precipitaciones, 4), data = datos)
summary(modelo_poli4)##
## Call:
## lm(formula = ProduccionCarne ~ poly(Precipitaciones, 4), data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61572 -37601 -20842 11009 186320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63364 10898 5.814 3.45e-06 ***
## poly(Precipitaciones, 4)1 53334 61646 0.865 0.395
## poly(Precipitaciones, 4)2 16606 61646 0.269 0.790
## poly(Precipitaciones, 4)3 -41319 61646 -0.670 0.508
## poly(Precipitaciones, 4)4 -75347 61646 -1.222 0.232
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 61650 on 27 degrees of freedom
## Multiple R-squared: 0.09287, Adjusted R-squared: -0.04152
## F-statistic: 0.6911 on 4 and 27 DF, p-value: 0.6046# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
# -----------------------------------------------------------------------------
limites <- range(datos$Precipitaciones)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(Precipitaciones = nuevos_puntos)
# PREDICCIÓN DE LA VARIABLE RESPUESTA Y DEL ERROR ESTÁNDAR
# -----------------------------------------------------------------------------
predicciones <- predict(modelo_poli4, newdata = nuevos_puntos, se.fit = TRUE,
level = 0.95)
# CÁLCULO DEL INTERVALO DE CONFIANZA SUPERIOR E INFERIOR 95%
# -----------------------------------------------------------------------------
intervalo_conf <- data.frame(inferior = predicciones$fit -
1.96*predicciones$se.fit,
superior = predicciones$fit +
1.96*predicciones$se.fit)
attach(datos)
plot(x =Precipitaciones, y = ProduccionCarne, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: Produccion de carne ~ Precipitación")
points(x = nuevos_puntos$Precipitaciones, predicciones$fit, col = "red", pch = 20)
points(x = nuevos_puntos$Precipitaciones, intervalo_conf$inferior, col = "blue", pch = 4)
points(x = nuevos_puntos$Precipitaciones, intervalo_conf$superior, col = "blue", pch = 4)
attach(datos)## The following objects are masked from datos (pos = 3):
##
## Entidad, PesoCarne, PrecioCarne, PrecioLeche, Precipitaciones,
## ProduccionCarne, ProduccionLeche, ValorCarne, ValorLecheplot(x = Precipitaciones, y = ProduccionCarne, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: Produccion de carne ~ Precipitación")
lines(x = nuevos_puntos$Precipitaciones, predicciones$fit, col = "red", lwd = 2)
lines(x = nuevos_puntos$Precipitaciones, intervalo_conf$inferior, col = "blue", lwd = 2)
lines(x = nuevos_puntos$Precipitaciones, intervalo_conf$superior, col = "blue", lwd = 2)
library(ggplot2)
ggplot(data = datos, aes(x = Precipitaciones, y = ProduccionCarne)) +
geom_point(color = "grey30", alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), color = "red") +
labs(title = "Polinomio de grado 4:Produccion de carne ~ Precipitación") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Comparación de modelos por contraste de hipótesis ANOVA
Identificar el modelo polinómico más simple que permite explicar la relación entre variables equivale a identificar el grado de polinomio a partir del cual ya no hay una mejora significativa del ajuste. Cuando se comparan dos modelos anidados (el modelo de menor tamaño está formado por un subset de predictores del modelo mayor), se puede saber si el modelo mayor aporta una mejora sustancial estudiando si los coeficientes de regresión de los predictores adicionales son distintos a cero. El test estadístico empleado para hacerlo es el ANOVA.
\[ Modelo_{menor}: \ \ y = \beta_0 + \beta_1x_1 +...+ \beta_kx_k \]
\[ Modelo_{mayor}: \ \ y = \beta_0 + \beta_1x_1 +...+ \beta_kx_k + \beta_{k+1}x_{k+1} + ... + \beta_{p}x_{p} \]
La hipótesis a contrastar es que todos los coeficientes de regresión de los predictores adicionales son igual a cero, frente a la hipótesis alternativa de que al menos uno es distinto.
\[ H_0: \beta_{k+1}= ... = \beta_{p} \]
El estadístico empleado es:
\[ F = \frac{(SEE_{Modelo_{menor}} - SEE_{Modelo_{mayor}})/(p-k)}{SEE_{Modelo_{mayor}}/(n-p-1)} \]
# Se ajustan modelos polinómicos de grado 1 a 5
# ----------------------------------------------------------------------------
modelo_1 <- lm(ProduccionCarne ~ Precipitaciones, data = datos)
modelo_2 <- lm(ProduccionCarne ~ poly(Precipitaciones, 2), data = datos)
modelo_3 <- lm(ProduccionCarne ~ poly(Precipitaciones, 3), data = datos)
modelo_4 <- lm(ProduccionCarne ~ poly(Precipitaciones, 4), data = datos)
modelo_5 <- lm(ProduccionCarne ~ poly(Precipitaciones, 5), data = datos)
anova(modelo_1, modelo_2, modelo_3, modelo_4, modelo_5)## Analysis of Variance Table
##
## Model 1: ProduccionCarne ~ Precipitaciones
## Model 2: ProduccionCarne ~ poly(Precipitaciones, 2)
## Model 3: ProduccionCarne ~ poly(Precipitaciones, 3)
## Model 4: ProduccionCarne ~ poly(Precipitaciones, 4)
## Model 5: ProduccionCarne ~ poly(Precipitaciones, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 1.1027e+11
## 2 29 1.0999e+11 1 275770301 0.0702 0.7931
## 3 28 1.0828e+11 1 1707235709 0.4347 0.5155
## 4 27 1.0261e+11 1 5677199202 1.4455 0.2401
## 5 26 1.0211e+11 1 494714631 0.1260 0.7255El p-value de la comparación entre el modelo lineal (modelo_1) y el cuadrático (modelo_2) es prácticamente cero, indicando que el modelo lineal no es suficiente. La comparación entre el modelo cuadrático y el cúbico también ha resultado en un p-value muy pequeño, evidencia suficiente de que el modelo cúbico es superior. En base a los resultados del ANOVA, el modelo cúbico es el mejor. Polinomios superiores no aportan una mejora significativa y polinomios inferiores pierden mucha capacidad de ajuste. Es importante recordar que las comparaciones por ANOVA pueden hacerse entre cualquier par de modelos, siempre y cuando estos sean anidados.
Comparación de modelos por Cross-Validation
Otra forma de identificar con que polinomio se consigue el mejor modelo es mediante cross-validation. El proceso consiste en ajustar un modelo para cada grado de polinomio y estimar su test error (Mean Square Error). El mejor modelo es aquel a partir del cual ya no hay una reducción sustancial del test error.
library(boot)
cv_MSE_k10 <- rep(NA,10)
for (i in 1:10) {
modelo <- glm(ProduccionCarne ~ poly(Precipitaciones, i), data = datos)
set.seed(17)
cv_MSE_k10[i] <- cv.glm(data = datos, glmfit = modelo, K = 10)$delta[1]
}## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 11.000000p4 <- ggplot(data = data.frame(polinomio = 1:10, cv_MSE = cv_MSE_k10),
aes(x = polinomio, y = cv_MSE)) +
geom_point(colour = c("firebrick3")) +
geom_path()
p4 <- p4 + theme(panel.grid.major = element_line(colour = 'gray90'))
p4 <- p4 + theme(plot.title = element_text(face = 'bold'))
p4 <- p4 + theme(panel.background = element_rect(fill = 'gray98'))
p4 <- p4 + labs(title = 'Test Error ~ Grado del polinomio')
p4 <- p4 + scale_x_continuous(breaks = 1:10)
p4
El método de cross-validation k=10 también indica que el mejor modelo es el que emplea un polinomio de grado 3.