Regresión Polinómica
En este apartado se realiza un análisis donde se efectúa la regresión no lineal para las variables PM10 (Concentración del aire) y espacio de supermercados y farmacias, donde en análisis anteriores se determinó que estas dos variables no tenían una correlación lineal.
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.18625datatable(aire)poly(aire$parques, degree = 4, raw = TRUE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 6 36 216 1296
## [2,] -4 16 -64 256
## [3,] -5 25 -125 625
## [4,] 2 4 8 16
## [5,] 8 64 512 4096poly(aire$parques, degree = 4, raw = FALSE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 0.2306840 0.147558771 0.10448042 0.02483518
## [2,] 0.1764430 0.003251611 -0.20663749 -0.17970681
## [3,] 0.1710189 -0.008192023 -0.22238206 -0.17410419
## [4,] 0.2089876 0.083318638 -0.05625387 -0.12786003
## [5,] 0.2415323 0.182937472 0.20535568 0.15007199Cálculo de Modelo Polinómico de Grado 4
modelo_poli4 <- lm(PM10 ~ poly(parques, 4), data = aire)
summary(modelo_poli4)##
## Call:
## lm(formula = PM10 ~ poly(parques, 4), data = aire)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.5174 -4.5737 -0.6405 5.3220 30.3574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.7529 0.8237 34.905 < 2e-16 ***
## poly(parques, 4)1 -13.3651 8.3195 -1.606 0.11142
## poly(parques, 4)2 -17.4246 8.3195 -2.094 0.03883 *
## poly(parques, 4)3 24.8361 8.3195 2.985 0.00359 **
## poly(parques, 4)4 6.7227 8.3195 0.808 0.42103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.319 on 97 degrees of freedom
## Multiple R-squared: 0.1456, Adjusted R-squared: 0.1104
## F-statistic: 4.133 on 4 and 97 DF, p-value: 0.003905# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
# -----------------------------------------------------------------------------
limites <- range(aire$parques)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(parques = 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 = parques, y = PM10, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: aire ~ parques")
points(x = nuevos_puntos$parques, predicciones$fit, col = "red", pch = 20)
points(x = nuevos_puntos$parques, intervalo_conf$inferior, col = "blue", pch = 4)
points(x = nuevos_puntos$parques, intervalo_conf$superior, col = "blue", pch = 4)
library(ggplot2)
ggplot(data = aire, aes(x = parques, 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: wage ~ age") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
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.
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
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
El método de cross-validation k=10 también indica que el mejor modelo es el que emplea un polinomio de grado 3.
Conclusión
Se determinó que el análisis para estas variables es del modelo 3, (modelo cúbico), cumpliendo con el p-value mayor de 0.05.