Tema 6: Modelos no lineales: Regresión Polinómica, Splines y GAMS
Diego Calderón
Mayo 2023
Introducción
En primer lugar cargamos la librería que contiene base de datos Clothing así como el resto de librerías que necesitaremos para resolver los ejercicios
library(Ecdat)## Loading required package: Ecfun##
## Attaching package: 'Ecfun'## The following object is masked from 'package:base':
##
## sign##
## Attaching package: 'Ecdat'## The following object is masked from 'package:datasets':
##
## Orangelibrary(splines)
library(gam)## Loading required package: foreach## Loaded gam 1.22library(akima)Para que sea más fácil trabajar con la base de datos, aplicamos la instrucción attach()
attach(Clothing)Ahora vamos a revisar los nombres de la base de datos y las características principales
names(Clothing)## [1] "tsales" "sales" "margin" "nown" "nfull" "npart" "naux"
## [8] "hoursw" "hourspw" "inv1" "inv2" "ssize" "start"summary(Clothing)## tsales sales margin nown
## Min. : 50000 Min. : 300 Min. :16.00 Min. : 1.000
## 1st Qu.: 495340 1st Qu.: 3904 1st Qu.:37.00 1st Qu.: 1.000
## Median : 694227 Median : 5279 Median :39.00 Median : 1.000
## Mean : 833584 Mean : 6335 Mean :38.77 Mean : 1.284
## 3rd Qu.: 976817 3rd Qu.: 7740 3rd Qu.:41.00 3rd Qu.: 1.295
## Max. :5000000 Max. :27000 Max. :66.00 Max. :10.000
## nfull npart naux hoursw
## Min. :1.000 Min. :1.000 Min. :1.000 Min. : 32.0
## 1st Qu.:1.923 1st Qu.:1.283 1st Qu.:1.333 1st Qu.: 80.0
## Median :1.956 Median :1.283 Median :1.367 Median :104.0
## Mean :2.069 Mean :1.566 Mean :1.390 Mean :121.1
## 3rd Qu.:2.066 3rd Qu.:2.000 3rd Qu.:1.367 3rd Qu.:145.2
## Max. :8.000 Max. :9.000 Max. :4.000 Max. :582.0
## hourspw inv1 inv2 ssize
## Min. : 5.708 Min. : 1000 Min. : 350 Min. : 16
## 1st Qu.:13.541 1st Qu.: 20000 1st Qu.: 10000 1st Qu.: 80
## Median :17.745 Median : 22207 Median : 22860 Median : 120
## Mean :18.955 Mean : 58257 Mean : 27829 Mean : 151
## 3rd Qu.:24.303 3rd Qu.: 62269 3rd Qu.: 22860 3rd Qu.: 190
## Max. :43.326 Max. :1500000 Max. :400000 Max. :1214
## start
## Min. :16.00
## 1st Qu.:37.00
## Median :40.00
## Mean :42.81
## 3rd Qu.:42.00
## Max. :90.00Regresión polinómica
Ejercicio: Una regresión polinómica de tercer grado que permita predecir el volumen de ventas anuales (tsales) en función de la inversión en locales comerciales (inv1).
Antes de comenzar el ejercicio, y para concer cual sería realmente el grado de polinomio a utilizar, veremos el test de hipotesis. Para ello ajustaremos todos los modleos que van hasta el polinomio de grado 5 y determinaremos el modelo más simple que sea sufiencite para explicar la relación entre tsales y inv. Para ello realizaremos un análisis de la varianza (ANOVA, utilizando un test F)
fit.1 <- lm(tsales ~ inv1, data = Clothing)
fit.2 <- lm(tsales ~ poly(inv1, 2), data = Clothing)
fit.3 <- lm(tsales ~ poly(inv1, 3), data = Clothing)
fit.4 <- lm(tsales ~ poly(inv1, 4), data = Clothing)
fit.5 <- lm(tsales ~ poly(inv1, 5), data = Clothing)
anova(fit.1, fit.2, fit.3, fit.4, fit.5)Como podemos observar, lo más eficiente sería usar un polinomio de grado 2, ya que el modelo cúbico 3 presenta un p-valor de 0.91.
En cualqueir caso, seguiremos las instrucciones del enunciado
Ajuste del modelo
En primer lugar, ajustamos el modelo utilizando la función lm utilizando el polinomio de tercer grado en inv1 (ploy(inv1,3).
fit <- lm(tsales ~ poly(inv1, 3), data = Clothing)
coef(summary(fit))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 833584.21 28379.71 29.3725407 1.738756e-101
## poly(inv1, 3)1 2220675.50 567594.22 3.9124350 1.075213e-04
## poly(inv1, 3)2 -1822025.35 567594.22 -3.2100844 1.435174e-03
## poly(inv1, 3)3 -61555.13 567594.22 -0.1084492 9.136943e-01Otras formas de ajustar el modelo
Otra forma de ajustar el modelo sería la siguiente:
fit2 <- lm(tsales ~ inv1 + I(inv1^2) + I(inv1^3), data = Clothing)
coef(fit2)## (Intercept) inv1 I(inv1^2) I(inv1^3)
## 7.289332e+05 2.032759e+00 -6.156011e-07 -4.196202e-13De una manera más compacta usando la función cbind()
fit2b <- lm(tsales ~ cbind(inv1, inv1^2, inv1^3), data = Clothing)
coef(fit2b)## (Intercept) cbind(inv1, inv1^2, inv1^3)inv1
## 7.289332e+05 2.032759e+00
## cbind(inv1, inv1^2, inv1^3) cbind(inv1, inv1^2, inv1^3)
## -6.156011e-07 -4.196202e-13Predicciones
A continuación creamos una cuadrícula de valores para poder hacer las predicciones. Usarmeos la función habitual de predict() especificando que queremos saber los errores estándar (se = TRUE)
agelims <- range(inv1)
inv1.grid <- seq(from = agelims[1], to = agelims[2])
preds <- predict(fit, newdata = list(inv1 = inv1.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)Finalmente visualizamos los datos y añadimos el ajuste polinómico de grado 3
plot(inv1, tsales, xlim = agelims, cex = .5, col = "darkgrey")
title("Regresión Polinómica de tercer grado")
lines(inv1.grid, preds$fit, lwd = 3, col = "blue")
matlines(inv1.grid, se.bands, lwd = 1, col = "red", lty = 3)
Tal como observaremos a continuación, que se produzca o no un conjunto ortogan de funciones no afectará al modelo de manera significativa ya que los valores obtenidos son casi idénticos (diferencia = 1.389068e-06)
preds2 <- predict(fit2, newdata = list(inv1 = inv1.grid), se = TRUE)
max(abs(preds$fit - preds2$fit))## [1] 1.389068e-06Modelo smooth spline
Ejercicio: Un modelo smooth spline que permita predecir el volumen de ventas anuales (tsales) en función de la inversión en automatización de procesos (inv2). Para ello, utilice un valor para λ=10, obtenga también el valor de λ por validación cruzada y represente gráficamente los resultados obtenidos.
En primer lugar, buscaremos los nudos ( Kknots ) que corresponden a los percentiles 25, 50 y 75 de la variable _inv_2
attr(bs(inv2 , df = 6), "knots")## 25% 50% 75%
## 10000.00 22859.85 22859.85Eejmplo con modelo spline
A modo de ejemplo, primero realizaremos un modelo spline
Para ajustar el modelo usamos la función bs. Como se observa en la gráfica, resulta poco natural, por lo que ya nos da una orientación de la necesidad de usar un modelo smooth spline
fit <- lm(tsales ~ bs(inv2, knots = c(10000.00, 22859.85, 22859.85 )), data = Clothing)
agelims <- range(inv2)
inv2.grid <- seq(from = agelims[1], to = agelims[2])
pred <- predict(fit, newdata = list(inv2 = inv2.grid), se = T)
plot(inv2, tsales, col = "gray")
title("modelo Spline")
lines(inv2.grid, pred$fit, lwd = 2, col = "blue")
lines(inv2.grid, pred$fit + 2 * pred$se, lty = "dashed", col = "red")
lines(inv2.grid, pred$fit - 2 * pred$se, lty = "dashed", col = "red")
Ajuste del modelo smooth spline
Para ajustar el modelo smooth spline usaremos la función smooth.splnie(). En este caso detemrinaremos que el valor λ conduzca a 10 grados de libertad ( df=10)
plot(inv2, tsales, xlim = agelims, cex = .5, col = "darkgrey")
title("Modelo Smooth Spline con λ = 10")
fit <- smooth.spline(inv2, tsales, df = 10)
lines(fit, col = "red", lwd = 2)
legend("topright", legend = c("10 DF"), col = c("red"), lty = 1, lwd = 2, cex = .8)
A continuación calculamos seleccionamos el nivel de suavidad por validación cruzada. Esto nos resulta en un valor de λ que produce 2.79 grados de libertad
fit2 <- smooth.spline(inv2, tsales, cv = TRUE)## Warning in smooth.spline(inv2, tsales, cv = TRUE): cross-validation with
## non-unique 'x' values seems doubtfulfit2$df## [1] 2.791762En el siguiente gráfico mostramos ambas opciones (λ = 10 en rojo y λ = 2.79 en azul)
plot(inv2, tsales, xlim = agelims, cex = .5, col = "darkgrey")
title("Modelo Smooth Spline con λ = 10 y λ = 2.79")
lines(fit, col = "red", lwd = 2)
lines(fit2, col = "blue", lwd = 2)
legend("topright", legend = c("10 DF", "2.79 DF"),
col = c("red", "blue"), lty = 1, lwd = 2, cex = .8)
Modelo GAM ( Generalized Additive Models )
Ejercicio: Un modelo GAM que permita predecir el volumen de ventas anuales (tsales) en función del número total de horas trabajadas (hoursw), el número de trabajadores a tiempo completo (nfull) y a tiempo parcial (npart), y la superficie destinada a ventas (ssize).
Ajuste del modelo lineal
En primer lugar ajustaremos el modelo pensandolo como un gram modelo de regresión lineal usando splines naturales. En este caso aplicaremos cuatro grados de libertad a cada parametro.
gam1 <- lm(tsales ~ ns(hoursw, 5) + ns(nfull, 5) +
ns(npart, 5) + ns(ssize, 5), data = Clothing)A continuación observamos el gráfico de las cuatro variables independientes propuestas. Como podemos observar, a pesar de ser una regresión lineal, la función plot.Gam() reconoce gam1 como de clse GAM a pesar de no serlo.
par(mfrow = c(1, 4))
plot.Gam(gam1 , se = TRUE , col = "red")
Ajuste del modelo con splines suavizados
A continuación ajustamos el modelo usando splines suavizados
gam.m3 <- gam(tsales ~ s(hoursw, 4) + s(nfull, 4) +
s(npart, 4) + s(ssize, 4), data = Clothing)Para verlo gráficamente usamos la función polt()
par(mfrow = c(1, 4))
plot(gam.m3, se = TRUE, col = "blue")
Justificación del uso del tip de GAM
Observamos que la variable hoursw puede ser lineal. Podemos realizar una serie de test ANOVA para para determinar cuál de estos tres modelos es el mejor: un GAM que excluye dicah variable hoursw (M1), un GAM que utiliza una función lineal de hoursw (M2), o un GAM que utiliza una función spline de hoursw (M3).
gam.m1 <- gam(tsales ~ s(nfull, 4) +
s(npart, 4) + s(ssize, 4), data = Clothing)
gam.m2 <- gam(tsales ~ hoursw + s(nfull, 4) +
s(npart, 4) + s(ssize, 4), data = Clothing)
anova(gam.m1, gam.m2, gam.m3, test = "F")Observamos que, si bien el modelo que utiliza una función lineal de hoursw (M2) es convincente, sí hay pruebas de que utilizar una función spline de hoursw (M3) es lo más adeucado.
Con la finción summary() podemos observar el ajuste de este GAM
summary(gam.m3)##
## Call: gam(formula = tsales ~ s(hoursw, 4) + s(nfull, 4) + s(npart,
## 4) + s(ssize, 4), data = Clothing)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -965434 -163583 -41494 141512 1929643
##
## (Dispersion Parameter for gaussian family taken to be 96028468879)
##
## Null Deviance: 1.358316e+14 on 399 degrees of freedom
## Residual Deviance: 3.677893e+13 on 383.0002 degrees of freedom
## AIC: 11268.94
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(hoursw, 4) 1 6.5460e+13 6.5460e+13 681.6724 < 2.2e-16 ***
## s(nfull, 4) 1 7.7740e+12 7.7740e+12 80.9551 < 2.2e-16 ***
## s(npart, 4) 1 1.2339e+12 1.2339e+12 12.8496 0.0003812 ***
## s(ssize, 4) 1 7.4996e+11 7.4996e+11 7.8098 0.0054577 **
## Residuals 383 3.6779e+13 9.6028e+10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(hoursw, 4) 3 7.8140 4.482e-05 ***
## s(nfull, 4) 3 22.0602 3.380e-13 ***
## s(npart, 4) 3 4.2314 0.005841 **
## s(ssize, 4) 3 11.7090 2.348e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los p-valores del “Anova para efectos paramétricos” (Anova for Parametric Effects) demuestran claramente que hoursw, nfull, npart y ssize son altamente significativos desde el punto de vista estadístico.
Otra forma d eobservarlo sería la siguiente:
gam.lo <- gam(tsales ~ s(hoursw, 4) + s(nfull, 4) +
s(npart, 4) + s(ssize, 4), data = Clothing)
par(mfrow = c(1,4))
plot(gam.lo, se = TRUE, col = "green")
Modelo con variables interactuando
Tal como hemos observado, las variables nfull y npart se comporetan de forma similar, por tanto podemos proponer un modelo en el cual interactuen entre ellas.
gam.lo.i <- gam(tsales ~ s(hoursw, 4) + lo(nfull, npart, span = 0.5) + s(ssize, 4), data = Clothing)
plot(gam.lo.i)## Warning in sqrt(out$se.fit[, TT]^2 + TS * object$var[, TT]): NaNs produced
summary(gam.lo.i)##
## Call: gam(formula = tsales ~ s(hoursw, 4) + lo(nfull, npart, span = 0.5) +
## s(ssize, 4), data = Clothing)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -975490 -170309 -47497 157311 1648554
##
## (Dispersion Parameter for gaussian family taken to be 98096439583)
##
## Null Deviance: 1.358316e+14 on 399 degrees of freedom
## Residual Deviance: 3.749936e+13 on 382.2703 degrees of freedom
## AIC: 11278.16
##
## Number of Local Scoring Iterations: 1
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(hoursw, 4) 1.00 5.9705e+13 5.9705e+13 608.637 < 2e-16 ***
## lo(nfull, npart, span = 0.5) 2.00 1.4438e+13 7.2188e+12 73.588 < 2e-16 ***
## s(ssize, 4) 1.00 3.6384e+11 3.6384e+11 3.709 0.05486 .
## Residuals 382.27 3.7499e+13 9.8096e+10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(hoursw, 4) 3.0 1.4917 0.2163
## lo(nfull, npart, span = 0.5) 6.7 16.4717 < 2.2e-16 ***
## s(ssize, 4) 3.0 13.5464 2.044e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusiones
En este capítulo hemos usado diversas formas de estimar una regresión cuando el modelo no responde a una forma lineal. A nivel personal he de decir que he aprovechado para usar GAM en unos datos de una investigación mía (https://www.data-in-brief.com/article/S2352-3409(23)00135-X/fulltext) y he podido probar su funcionalidad, ya que al realizar en un primer momento una regresión líneal no obtuve relaciones significativas y al usar Splines suavizados sí