GAM REGRESSION
2023-10-12
#1.- INTRODUCCIÓN
Para entender mejor el contenido de este capítulo es recomendable leer primero: Introducción a la Regresión Lineal Múltiple, Ejemplo práctico de regresión lineal simple, múltiple, polinomial e interacción entre predictores y Regresión logística simple y múltiple.
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.
Regression splines: Se trata de una extensión de la regresión polinómica y de las step functions que consigue una mayor flexibilidad. Consiste en dividir el rango del predictor X en K subintervalos. Para cada una de las nuevas regiones se ajusta una función polinómica, introduciendo una serie de restricciones que hacen que los extremos de cada función se aproximen a los de las funciones de las regiones colindantes.
Smoothing splines: El concepto es similar a regression splines pero consigue la aproximación de los extremos de las funciones colindantes de forma distinta. Local regression: Se asemeja a regression splines y smoothing splines en cuanto a que también se realizan ajustes por regiones, pero en este método las regiones solapan las unas con las otras.
Local regression: Se asemeja a regression splines y smoothing splines en cuanto a que también se realizan ajustes por regiones, pero en este método las regiones solapan las unas con las otras.
Generalized additive models: Es el resultado de extender los métodos anteriores para emplear múltiples predictores.
##2.- EJERCICIO
El set de datos Wage del paquete ISRL contiene información sobre 3000 trabajadores. Entre las 12 variables registradas se encuentra el salario (wage) y la edad (age). Dada la relación no lineal existente entre estas dos variables, se recurre a un modelo polinómico de grado 4 que permita predecir el salario en función de la edad.
La función lm() permite ajustar modelos lineales por mínimos cuadrados. En el argumento formula se especifica la variable dependiente y los predictores, que en este caso son wage ~ age + age^2 + age^3 + age^4. lm() tiene flexibilidad a la hora de interpretar el argumento fórmula, por lo que existen diferentes formas de obtener el mismo ajuste. Todas las que se muestran a continuación son equivalentes:
lm(wage ~ age + I(age^2) + I(age^3) + I(age^4), data = Wage)
lm(wage ~ cbind(age + age^2 + age^3 + age^4), data = Wage)
lm(wage ~ poly(age, 4, raw = TRUE), data = Wage)
lm(wage ~ poly(age, 4), data = Wage)
library(ISLR)## Warning: package 'ISLR' was built under R version 4.2.3data("Wage")poly(Wage$age, degree = 4, raw = T, simple = T)[1:5,]## 1 2 3 4
## [1,] 18 324 5832 104976
## [2,] 24 576 13824 331776
## [3,] 45 2025 91125 4100625
## [4,] 43 1849 79507 3418801
## [5,] 50 2500 125000 6250000poly(Wage$age, degree = 4, raw = F, simple = T)[1:5,]## 1 2 3 4
## [1,] -0.0386247992 0.055908727 -0.0717405794 0.08672985
## [2,] -0.0291326034 0.026298066 -0.0145499511 -0.00259928
## [3,] 0.0040900817 -0.014506548 -0.0001331835 0.01448009
## [4,] 0.0009260164 -0.014831404 0.0045136682 0.01265751
## [5,] 0.0120002448 -0.009815846 -0.0111366263 0.01021146# Modelo polinómico
# ========================
library(polynom)## Warning: package 'polynom' was built under R version 4.2.3modelo_poli4 <- lm(wage ~ poly(age,4), data = Wage)
summary(modelo_poli4)##
## Call:
## lm(formula = wage ~ poly(age, 4), data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -98.707 -24.626 -4.993 15.217 203.693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.7036 0.7287 153.283 < 2e-16 ***
## poly(age, 4)1 447.0679 39.9148 11.201 < 2e-16 ***
## poly(age, 4)2 -478.3158 39.9148 -11.983 < 2e-16 ***
## poly(age, 4)3 125.5217 39.9148 3.145 0.00168 **
## poly(age, 4)4 -77.9112 39.9148 -1.952 0.05104 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.91 on 2995 degrees of freedom
## Multiple R-squared: 0.08626, Adjusted R-squared: 0.08504
## F-statistic: 70.69 on 4 and 2995 DF, p-value: < 2.2e-16library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(modelo_poli4, type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## wage
## -----------------------------------------------
## poly(age, 4)1 447.068***
## (39.915)
##
## poly(age, 4)2 -478.316***
## (39.915)
##
## poly(age, 4)3 125.522***
## (39.915)
##
## poly(age, 4)4 -77.911*
## (39.915)
##
## Constant 111.704***
## (0.729)
##
## -----------------------------------------------
## Observations 3,000
## R2 0.086
## Adjusted R2 0.085
## Residual Std. Error 39.915 (df = 2995)
## F Statistic 70.689*** (df = 4; 2995)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3ggplot(data = Wage, aes(x = age, y = wage)) +
geom_point(color = "grey30", alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), color = "red") +
labs(title = "Polinomio de grado 4: wage ~ age") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Cuando se realiza regresión polinómica, se debe decidir qué grado de
polinomio emplear. Cuanto mayor sea el polinomio más flexibilidad tendrá
el modelo pero, a su vez, más riesgo de overfitting. Acorde al principio
de parsimonia, el grado óptimo es el grado más bajo que permita explicar
la relación entre ambas variables. Para identificarlo se puede recurrir
a dos estrategias distintas: contraste de hipótesis o
cross-validation.
Comparación de los modelos
m1 <- lm(wage ~age, data = Wage)
m2 <- lm(wage ~poly(age, 2), data = Wage)
m3 <- lm(wage ~poly(age, 3), data = Wage)
m4 <- lm(wage ~poly(age, 4), data = Wage)
m5 <- lm(wage ~poly(age, 5), data = Wage)
m6 <- lm(wage ~poly(age, 6), data = Wage)library(stargazer)stargazer(m1, m2, m3,m4,m5,m6, type = "text")##
## ===========================================================================================================================================================================
## Dependent variable:
## -------------------------------------------------------------------------------------------------------------------------------------------------------
## wage
## (1) (2) (3) (4) (5) (6)
## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## age 0.707***
## (0.065)
##
## poly(age, 2)1 447.068***
## (39.993)
##
## poly(age, 2)2 -478.316***
## (39.993)
##
## poly(age, 3)1 447.068***
## (39.933)
##
## poly(age, 3)2 -478.316***
## (39.933)
##
## poly(age, 3)3 125.522***
## (39.933)
##
## poly(age, 4)1 447.068***
## (39.915)
##
## poly(age, 4)2 -478.316***
## (39.915)
##
## poly(age, 4)3 125.522***
## (39.915)
##
## poly(age, 4)4 -77.911*
## (39.915)
##
## poly(age, 5)1 447.068***
## (39.916)
##
## poly(age, 5)2 -478.316***
## (39.916)
##
## poly(age, 5)3 125.522***
## (39.916)
##
## poly(age, 5)4 -77.911*
## (39.916)
##
## poly(age, 5)5 -35.813
## (39.916)
##
## poly(age, 6)1 447.068***
## (39.906)
##
## poly(age, 6)2 -478.316***
## (39.906)
##
## poly(age, 6)3 125.522***
## (39.906)
##
## poly(age, 6)4 -77.911*
## (39.906)
##
## poly(age, 6)5 -35.813
## (39.906)
##
## poly(age, 6)6 62.708
## (39.906)
##
## Constant 81.705*** 111.704*** 111.704*** 111.704*** 111.704*** 111.704***
## (2.846) (0.730) (0.729) (0.729) (0.729) (0.729)
##
## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## Observations 3,000 3,000 3,000 3,000 3,000 3,000
## R2 0.038 0.082 0.085 0.086 0.087 0.087
## Adjusted R2 0.038 0.081 0.084 0.085 0.085 0.085
## Residual Std. Error 40.929 (df = 2998) 39.993 (df = 2997) 39.933 (df = 2996) 39.915 (df = 2995) 39.916 (df = 2994) 39.906 (df = 2993)
## F Statistic 119.312*** (df = 1; 2998) 134.004*** (df = 2; 2997) 92.894*** (df = 3; 2996) 70.689*** (df = 4; 2995) 56.708*** (df = 5; 2994) 47.692*** (df = 6; 2993)
## ===========================================================================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01anova(m1,m2,m3,m4, m5, m6)## Analysis of Variance Table
##
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Model 6: wage ~ poly(age, 6)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.6636 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8936 0.001675 **
## 4 2995 4771604 1 6070 3.8117 0.050989 .
## 5 2994 4770322 1 1283 0.8054 0.369565
## 6 2993 4766389 1 3932 2.4692 0.116201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aquí vemos que hay diferencias significativas.
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. Para una descripción detallada de cross-validation y del código empleado a continuación ver capítulo Validación de modelos de regresión: Cross-validation, OneLeaveOut, Bootstrap.
library(boot)library(boot)CV_MSE_k10 <- rep(NA,10)
for (i in 1:10) {
modelo <- glm(wage ~ poly(age, i), data = Wage)
set.seed(17)
CV_MSE_k10[i] <- cv.glm(data = Wage, glmfit = modelo, K = 10)$delta[1]
}
p4 <- 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