Regresión polinómica
Descarga de datos
- Para la descarga de este código
Para fines de reproducibilidad inmediata se incluye todo el código para su descarga
xfun::embed_file("U2A5.Rmd")- Para la descarga de datos utilizados en este codigo
Para fines de reproducibilidad inmediata se incluye todos los datos para su descarga
xfun::embed_file("Aire.xlsx")Librerías y Paquetes
library(pacman) #Para importar la biblioteca "pacman"
p_load("base64enc", "htmltools", "mime", "xfun", "prettydoc","readr", "knitr","DT","dplyr", "ggplot2","plotly", "gganimate","gifski","scales")Importar Datos
library(DT)
library(readxl)
aire <- read_excel("Aire.xlsx")
data("aire")## Warning in data("aire"): data set 'aire' not foundaire$PM10[1:5]## [1] 32.31750 28.50333 12.30250 26.29667 22.18625poly(aire$SupermercadosyFarmacias, degree = 4, raw = TRUE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 8 64 512 4096
## [2,] 4 16 64 256
## [3,] 0 0 0 0
## [4,] 4 16 64 256
## [5,] 4 16 64 256poly(aire$SupermercadosyFarmacias, degree = 4, raw = FALSE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 0.2043764 0.1160876 -0.05477833 -0.17704130
## [2,] 0.1705073 0.0387702 -0.14344878 -0.16021328
## [3,] 0.1366383 -0.0205399 -0.17030312 -0.07490379
## [4,] 0.1705073 0.0387702 -0.14344878 -0.16021328
## [5,] 0.1705073 0.0387702 -0.14344878 -0.16021328CÁLCULO DEL MODELO POLINÓMICO DE GRADO 4
modelo_poli4 <- lm(PM10 ~ poly(SupermercadosyFarmacias, 4), data = aire)
summary(modelo_poli4)##
## Call:
## lm(formula = PM10 ~ poly(SupermercadosyFarmacias, 4), data = aire)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.8069 -4.5032 -0.5977 5.5051 31.2660
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.7529 0.8436 34.086 <2e-16 ***
## poly(SupermercadosyFarmacias, 4)1 -11.1770 8.5195 -1.312 0.1926
## poly(SupermercadosyFarmacias, 4)2 -18.0481 8.5195 -2.118 0.0367 *
## poly(SupermercadosyFarmacias, 4)3 18.9606 8.5195 2.226 0.0284 *
## poly(SupermercadosyFarmacias, 4)4 2.7231 8.5195 0.320 0.7499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.519 on 97 degrees of freedom
## Multiple R-squared: 0.104, Adjusted R-squared: 0.0671
## F-statistic: 2.816 on 4 and 97 DF, p-value: 0.02933El p-value obtenido es bajo (0.02933), lo que indica que al menos uno de los predictores introducidos en el modelo está relacionado con la variable PM10. Los p-values individuales de cada predictor son todos muy bajos a excepción de Supermercados y farmacias^3 y elevado a la cuarta, lo que apunta a que un polinomio de grado 3 es suficiente para modelar el salario en función de la edad.
Si se quiere generar una representación gráfica del modelo ajustado:
\[ PM10=28.7529−11.1770Supermercadosyfarmacias−18.0481Supermercadosyfarmacias^2+18.9606Supermercadosyfarmacias^3+2.7231Supermercadosyfarmacias^4 \]
# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
# -----------------------------------------------------------------------------
limites <- range(aire$SupermercadosyFarmacias)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(SupermercadosyFarmacias = 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(aire)
plot(x = SupermercadosyFarmacias, y = PM10, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: aire ~ SupermercadosyFarmacias")
points(x = nuevos_puntos$SupermercadosyFarmacias, predicciones$fit, col = "red", pch = 20)
points(x = nuevos_puntos$SupermercadosyFarmacias, intervalo_conf$inferior, col = "blue", pch = 4)
points(x = nuevos_puntos$SupermercadosyFarmacias, intervalo_conf$superior, col = "blue", pch = 4)
attach(aire)## The following objects are masked from aire (pos = 3):
##
## EstacionesdeTransporte, LugaresdeTrabajo, NO2, O3, parques, PM10,
## SO2, SupermercadosyFarmacias, TiendasyOcio, ZonasResidencialesplot(x = SupermercadosyFarmacias, y = PM10, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: PM10 ~ SupermercadosyFarmacias")
lines(x = nuevos_puntos$SupermercadosyFarmacias, predicciones$fit, col = "red", lwd = 2)
lines(x = nuevos_puntos$SupermercadosyFarmacias, intervalo_conf$inferior, col = "blue", lwd = 2)
lines(x = nuevos_puntos$SupermercadosyFarmacias, intervalo_conf$superior, col = "blue", lwd = 2)
library(ggplot2)
ggplot(data = aire, aes(x = SupermercadosyFarmacias, y = PM10)) +
geom_point(color = "grey30", alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), color = "red") +
labs(title = "Polinomio de grado 4: Parques ~ PM10") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Se puede explicar que las concentraciones de PM10 aumentan cuando el índice de supermercados y farmacias se encientran entre el -20 y -10
COMPARACIÓN DE MODELOS POR CONTRASTE DE HIPÓTESIS ANOVA
# Se ajustan modelos polinómicos de grado 1 a 5
# ----------------------------------------------------------------------------
modelo_1 <- lm(PM10 ~ SupermercadosyFarmacias, data = aire)
modelo_2 <- lm(PM10 ~ poly(SupermercadosyFarmacias, 2), data = aire)
modelo_3 <- lm(PM10 ~ poly(SupermercadosyFarmacias, 3), data = aire)
modelo_4 <- lm(PM10 ~ poly(SupermercadosyFarmacias, 4), data = aire)
modelo_5 <- lm(PM10 ~ poly(SupermercadosyFarmacias, 5), data = aire)
anova(modelo_1, modelo_2, modelo_3, modelo_4, modelo_5)## Analysis of Variance Table
##
## Model 1: PM10 ~ SupermercadosyFarmacias
## Model 2: PM10 ~ poly(SupermercadosyFarmacias, 2)
## Model 3: PM10 ~ poly(SupermercadosyFarmacias, 3)
## Model 4: PM10 ~ poly(SupermercadosyFarmacias, 4)
## Model 5: PM10 ~ poly(SupermercadosyFarmacias, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 100 7733.0
## 2 99 7407.3 1 325.73 4.6822 0.03296 *
## 3 98 7047.8 1 359.50 5.1676 0.02524 *
## 4 97 7040.4 1 7.42 0.1066 0.74477
## 5 96 6678.5 1 361.84 5.2012 0.02478 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En base a los resultados del ANOVA, el modelo cúbico es el mejor ya que presenta p-value ligeramente grande (0.75). Polinomios superiores no aportan una mejora significativa y polinomios inferiores pierden mucha capacidad de ajuste.
COMPARACIÓN DE MODELOS POR CROSS-VALIDATION
library(boot)
cv_MSE_k10 <- rep(NA,10)
for (i in 1:10) {
modelo <- glm(PM10 ~ poly(SupermercadosyFarmacias, i), data = aire)
set.seed(17)
cv_MSE_k10[i] <- cv.glm(data = aire, 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
Conclusión
Gracias a esta actividad se concluyó que el modelo es del número tres debido a que eso se determinó gracias al análisis de las variables, además de que el p-value resulto ser mayor al mínimo ya estimado de 0.05