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.
- Descarga de este archivo
xfun::embed_file("U2A5.Rmd")library(DT)
library(readxl)
aire <- read_excel("Aire.xlsx")
datatable(aire)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. Acorde al R2, el modelo es capaz de explicar el 6.71% de la variabilidad observada en PM10 (es un % muy bajo).
Si se quiere generar una representación gráfica del modelo ajustado:
\[ PM10= 28.7529-11.1770Supermercados y farmacias-18.0481Supermercados y farmacias^2+18.9606Supermercados y farmacias^3 + 2.7231Supermercados y farmacias^4 \]
se tiene que representar la curva que describe su ecuación. Para conseguirlo, es necesario predecir el valor de la variable respuesta en puntos interpolados dentro del rango del predictor.
# 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.
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.
Se realizó un análisis de regresión no polinómica a partir de dos variables que no tenían mucha relación en pasados estudios.