Analisis de Regresión no Lineal
El analisis de regresión no lineal se fundamente de establecer una relación entre dos variables, en este caso de las colmenas y precipitacion nacional en México, con una aparente no linealidad entre las variables de estudio, donde dado a estudios se estableceria que no hay una correlación directa entre ambas.
setwd("~/estadistica aplicada")
library(DT)
library(pacman)
library(psych)
p_load("MASS", "ggplot2","readr", "prettydoc")
datos <- read_csv("abejas.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## TonPlagui = col_double(),
## Colmenas = col_double(),
## Deforestacion = col_double(),
## Tempmedia = col_double(),
## Densidadpoblacional = col_double()
## )datos$Colmenas[1:5]## [1] 1862372 1783854 1727234 1745078 1732112poly(datos$Deforestacion, degree = 4, raw = TRUE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] 79672 6347627584 5.057282e+14 4.029238e+19
## [2,] 191071 36508127041 6.975644e+15 1.332843e+21
## [3,] 185741 34499719081 6.408012e+15 1.190231e+21
## [4,] 135953 18483218209 2.512849e+15 3.416294e+20
## [5,] 170421 29043317241 4.949591e+15 8.435143e+20poly(datos$Deforestacion, degree = 4, raw = FALSE, simple = TRUE)[1:5,]## 1 2 3 4
## [1,] -0.36331411 0.4291656 -0.3667906 0.414947176
## [2,] -0.05762284 -0.2502034 0.1000163 0.202418750
## [3,] -0.07224895 -0.2387509 0.1315644 0.157600218
## [4,] -0.20887277 -0.0295821 0.2463577 -0.345672763
## [5,] -0.11428874 -0.1940524 0.2072683 0.002248779Cálculo del modelo polinómico de grado 4
modelo_polin4 <- lm(Colmenas ~ poly(Deforestacion, 4), data = datos)
summary(modelo_polin4)##
## Call:
## lm(formula = Colmenas ~ poly(Deforestacion, 4), data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -118694 -53241 -40602 32441 360580
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1850914 29559 62.618 <2e-16 ***
## poly(Deforestacion, 4)1 155046 125408 1.236 0.238
## poly(Deforestacion, 4)2 36969 125408 0.295 0.773
## poly(Deforestacion, 4)3 -52735 125408 -0.421 0.681
## poly(Deforestacion, 4)4 36112 125408 0.288 0.778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 125400 on 13 degrees of freedom
## Multiple R-squared: 0.1261, Adjusted R-squared: -0.1428
## F-statistic: 0.4688 on 4 and 13 DF, p-value: 0.7578Se obtuvo un valor de p-value de 0.7578 con referencia a las variables introducidas en el analisis.
Si se quiere generar una representación gráfica del modelo ajustado:
\[ Colmenas= 1850914+155046 Deforestacion+36969Deforestacion^2-52735Deforestacion^3 + 36112Deforestacion^4 \]
Con la realización de la prediccion de la valor de las varianles de respuesta en puntos interpolados dentro del rango predictor nos permitirá lograr una representacion de la curva que describa a la ecuacion identificada.
# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
# -----------------------------------------------------------------------------
limites <- range(datos$Deforestacion)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(Deforestacion = nuevos_puntos)
# PREDICCIÓN DE LA VARIABLE RESPUESTA Y DEL ERROR ESTÁNDAR
# -----------------------------------------------------------------------------
predicciones <- predict(modelo_polin4, newdata = nuevos_puntos, se.fit = TRUE,
level = 0.90)
# 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(datos)
plot(x = Deforestacion, y = Colmenas, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: Colmenas ~ Deforestacion")
points(x = nuevos_puntos$Deforestacion, predicciones$fit, col = "red", pch = 20)
points(x = nuevos_puntos$Deforestacion, intervalo_conf$inferior, col = "blue", pch = 4)
points(x = nuevos_puntos$Deforestacion, intervalo_conf$superior, col = "blue", pch = 4)
attach(datos)## The following objects are masked from datos (pos = 3):
##
## Colmenas, Deforestacion, Densidadpoblacional, Tempmedia, TonPlaguiplot(x = Deforestacion, y = Colmenas, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: Colmenas ~ Deforestacion")
lines(x = nuevos_puntos$Deforestacion, predicciones$fit, col = "red", lwd = 2)
lines(x = nuevos_puntos$Deforestacion, intervalo_conf$inferior, col = "blue", lwd = 2)
lines(x = nuevos_puntos$Deforestacion, intervalo_conf$superior, col = "blue", lwd = 2)
library(ggplot2)
ggplot(data = datos, aes(x = Deforestacion, y = Colmenas)) +
geom_point(color = "grey30", alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), color = "red") +
labs(title = "Polinomio de grado 4: Colmenas ~ Deforestacion") +
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 modelos por contraste de hipótesis ANOVA
# Se ajustan modelos polinómicos de grado 1 a 5
# ----------------------------------------------------------------------------
modelo_1 <- lm(Colmenas ~ Deforestacion, data = datos)
modelo_2 <- lm(Colmenas ~ poly(Deforestacion, 2), data = datos)
modelo_3 <- lm(Colmenas ~ poly(Deforestacion, 3), data = datos)
modelo_4 <- lm(Colmenas ~ poly(Deforestacion, 4), data = datos)
modelo_5 <- lm(Colmenas ~ poly(Deforestacion, 5), data = datos)
anova(modelo_1, modelo_2, modelo_3, modelo_4, modelo_5)## Analysis of Variance Table
##
## Model 1: Colmenas ~ Deforestacion
## Model 2: Colmenas ~ poly(Deforestacion, 2)
## Model 3: Colmenas ~ poly(Deforestacion, 3)
## Model 4: Colmenas ~ poly(Deforestacion, 4)
## Model 5: Colmenas ~ poly(Deforestacion, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 16 2.0991e+11
## 2 15 2.0854e+11 1 1366701693 0.0815 0.7802
## 3 14 2.0576e+11 1 2780931851 0.1658 0.6910
## 4 13 2.0445e+11 1 1304048254 0.0778 0.7851
## 5 12 2.0125e+11 1 3199969247 0.1908 0.6700A base de 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)##
## Attaching package: 'boot'## The following object is masked from 'package:psych':
##
## logitcv_MSE_k10 <- rep(NA,10)
for (i in 1:10) {
modelo <- glm(Colmenas ~ poly(Deforestacion, i), data = datos)
set.seed(17)
cv_MSE_k10[i] <- cv.glm(data = datos, glmfit = modelo, K = 10)$delta[1]
}## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000
## Warning in cv.glm(data = datos, glmfit = modelo, K = 10): 'K' has been set to
## 9.000000p4 <- 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
En el analisis de la regresion no lineal se busca la determinación de la relacion entre variables que aparentemente no cuentan con una correlación directa con respecto a la tendencia de sus datos registrados donde uno de los principales factores a tomar en cuenta son los valores obtenidos a partir de p-value. En este caso los datos a analisar se dieron entre las colmenas y la deforestacion donde 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.
Descargar Codigo&Datos
xfun::embed_file("u2a5-.rmd")xfun::embed_file("abejas.csv")