U2A4: Regresión lineal multiple.
Introducción
Para el siguiente trabajo hablaremos sobre la regresión lineal múltiple y estaremos tomando en cuenta como es que la papa actúa en los sembradíos cuando se siembra, cosecha, cando sucede siniestréo y el valor de la producción total. Obtener la regresión lineal múltiple puede ser de información bastante valiosa ya que nos indicará el como se relacionan unas cantidades con otras y poder llegar a una correlación mas amplia. Además de esto se puede diferir con la información que estamos a punto de procesar, causalidades o eventos futuros por los cuales guiarse en el futuro.
Las papas como las conocemos son un tubérculo que pertenece en varias comidas mexicanas y platillos caseros que se consumen por todo el país. Este tubérculo es muy bien conocido entre los agricultores por ocupar muy pocos recursos para abastecer de comida a la población en comparación a los demás. La papa se cultiva en 22 países de la república y es considerado en México como el 7mo cultivo más importante.
Valor
Importación de paquetes
library(pacman)
# Si esta en linux descarga unas bibliotecas necesarias para correr este código
if(Sys.info()[['sysname']]=='Linux'){
system('sudo -i apt-get install -y r-cran-psych r-cran-car r-cran-lmtest',input = rstudioapi::askForPassword("sudo password"))
}
p_load(GGally, ggplot2,readr,DT,rmdformats, xfun,data.table,psych, gridExtra, lmtest, corrplot, car,dplyr)Importacion de datos
names = c("AÑO","ENE","FEB","MAR","ABR","MAY","JUN","JUL","AGO","SEP","OCT","NOV","DIC","ACUMULADO","MEDIA","MESES")
lluvias <- read.table("lluvias.txt")
names(lluvias) <- names
temperaturas = read.table("temperaturas.txt")
names(temperaturas) <- names
papas <- read.csv("papasV2.csv")
papas <- papas[-c(seq(26,30)),]
papas$lluvias <- lluvias$MEDIA
papas$temperaturas <- temperaturas$MEDIADescarga de los datos
xfun::embed_file("papasV2.csv")xfun::embed_file("temperaturas.txt")xfun::embed_file("lluvias.txt")Primera visualizacion de los datos
datatable(papas)Visualización de tabla de predictores
papas <- data.table(papas[-1])
round(cor(x=papas, method="pearson"),3)## Superficie_Sembrada Superficie_Cosechada
## Superficie_Sembrada 1.000 0.937
## Superficie_Cosechada 0.937 1.000
## Superficie_Siniestrada 0.366 0.017
## Valor_de_la_Produccion 0.875 0.897
## lluvias -0.183 -0.241
## temperaturas 0.135 0.094
## Superficie_Siniestrada Valor_de_la_Produccion lluvias
## Superficie_Sembrada 0.366 0.875 -0.183
## Superficie_Cosechada 0.017 0.897 -0.241
## Superficie_Siniestrada 1.000 0.116 0.117
## Valor_de_la_Produccion 0.116 1.000 -0.379
## lluvias 0.117 -0.379 1.000
## temperaturas 0.136 -0.074 0.008
## temperaturas
## Superficie_Sembrada 0.135
## Superficie_Cosechada 0.094
## Superficie_Siniestrada 0.136
## Valor_de_la_Produccion -0.074
## lluvias 0.008
## temperaturas 1.000En base a los resultados de las correlaciones de pearson se podría decir que el valor de producción de la papa esta relacionada de una gran manera, aproximadamente un 95% con la superficie cosechada de la misma. y se puede ver que la relación que existe entre la superficie cosechada tiene una relación negativa con la superficie siniestrada, esto básicamente significa que si uno sube el otro baja, lo cual es lógico, porque la superficie siniestrada es todo el material que se pierde.
Visualización del comportamiento de los datos en histogramas de frecuencia.
#par (mfrow = c(2,2))
#hist(papas$Superficie_Sembrada)
#hist(papas$Superficie_Cosechada)
#hist(papas$Superficie_Siniestrada)
#hist(papas$Valor_de_la_Produccion)
multi.hist(x = papas, dcol = c("blue", "red"), dlty = c("dotted", "solid"), main = "")
Como podemos ver en la siguiente gráfica múltiple, se muestran que el valor de la producción y la superficie siniestrada no están correlacionadas tanto como la superficie sembrada y la superficie cosechada, esto puede ser un poco obvio al hablar de teoría pero esto nos ayuda a reafirmar la posición que se tiene sobre estas variables, además de que se pueden ver unos picos interesantes que hacen falta investigar para saber que fue lo que causo que hubiera tanta superficie siniestrada y tanto valor de producción en las primeras etapas de la gráfica
Matriz de correlación multiple
ggpairs(papas, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")
En está gráfica se puede apreciar como algunas de las gráficas están en efecto correlacionadas con otras, la manera de saber como es que están verdaderamente correlacionadas es cuando estas pasan el nivel de 0.05 en sus niveles, por lo que se puede apreciar que en almenos 3 tablas la correlación fue exitosa y alta, mientras que en las demás se puede ver que no tienen correlación en absoluto
Generación de modelo para el valor de producción
modelo <- lm(Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada + Superficie_Siniestrada + lluvias + temperaturas, data = papas )
summary(modelo)##
## Call:
## lm(formula = Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada +
## Superficie_Siniestrada + lluvias + temperaturas, data = papas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -661087 -202627 15553 155161 1342489
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4787216.42 2387853.91 2.005 0.0587 .
## Superficie_Sembrada 119.44 65.23 1.831 0.0820 .
## Superficie_Cosechada 146.54 70.66 2.074 0.0512 .
## Superficie_Siniestrada NA NA NA NA
## lluvias -31919.96 14007.98 -2.279 0.0338 *
## temperaturas -160143.93 73138.51 -2.190 0.0406 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 415200 on 20 degrees of freedom
## Multiple R-squared: 0.8766, Adjusted R-squared: 0.8519
## F-statistic: 35.51 on 4 and 20 DF, p-value: 8.008e-09Seleccion de los mejores predictores
En este caso se usará la estrategia stepwise mixto. El valor matematico empleado para determinar la calidad del modelo es Akaike(AIC)
step(object=modelo, direction="both", trace=1)## Start: AIC=651.25
## Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada +
## Superficie_Siniestrada + lluvias + temperaturas
##
##
## Step: AIC=651.25
## Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada +
## lluvias + temperaturas
##
## Df Sum of Sq RSS AIC
## <none> 3.4478e+12 651.25
## - Superficie_Sembrada 1 5.7792e+11 4.0258e+12 653.12
## - Superficie_Cosechada 1 7.4155e+11 4.1894e+12 654.12
## - temperaturas 1 8.2650e+11 4.2743e+12 654.62
## - lluvias 1 8.9514e+11 4.3430e+12 655.02##
## Call:
## lm(formula = Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada +
## lluvias + temperaturas, data = papas)
##
## Coefficients:
## (Intercept) Superficie_Sembrada Superficie_Cosechada
## 4787216.4 119.4 146.5
## lluvias temperaturas
## -31920.0 -160143.9Como se puede ver el mejor modelo ha sido la superficie cosechada.
modelo <- (lm(Valor_de_la_Produccion ~ Superficie_Cosechada + Superficie_Sembrada, data = papas))
summary(modelo)##
## Call:
## lm(formula = Valor_de_la_Produccion ~ Superficie_Cosechada +
## Superficie_Sembrada, data = papas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -861805 -320832 -37316 262615 1280565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -987406.46 234825.12 -4.205 0.000366 ***
## Superficie_Cosechada 193.10 80.68 2.393 0.025659 *
## Superficie_Sembrada 81.85 75.09 1.090 0.287462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 486000 on 22 degrees of freedom
## Multiple R-squared: 0.814, Adjusted R-squared: 0.7971
## F-statistic: 48.14 on 2 and 22 DF, p-value: 9.211e-09Se mostrará el intervalo de confianza para cada uno de los coeficientes parciales de la regresión, que en este caso solo es la superficie cosechada.
modelo <- lm(Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada + Superficie_Siniestrada + lluvias + temperaturas, data = papas )
confint(modelo)## 2.5 % 97.5 %
## (Intercept) -1.937596e+05 9768192.3977
## Superficie_Sembrada -1.663551e+01 255.5127
## Superficie_Cosechada -8.443367e-01 293.9314
## Superficie_Siniestrada NA NA
## lluvias -6.114008e+04 -2699.8319
## temperaturas -3.127082e+05 -7579.6755Validación de condiciones para la regresión multiple lineal
plot1 <- ggplot(papas, aes(Superficie_Sembrada, modelo$residuals))+
geom_point()+ geom_smooth(color="darkred")+geom_hline(yintercept = 0)+
theme_bw()
plot2 <- ggplot(papas, aes(Superficie_Cosechada, modelo$residuals))+
geom_point()+ geom_smooth(color="darkred")+geom_hline(yintercept = 0)+
theme_bw()
plot3 <- ggplot(papas, aes(Superficie_Siniestrada, modelo$residuals))+
geom_point()+ geom_smooth(color="darkred")+geom_hline(yintercept = 0)+
theme_bw()
plot4 <- ggplot(papas, aes(lluvias, modelo$residuals))+
geom_point()+ geom_smooth(color="darkred")+geom_hline(yintercept = 0)+
theme_bw()
plot5 <- ggplot(papas, aes(temperaturas, modelo$residuals))+
geom_point()+ geom_smooth(color="darkred")+geom_hline(yintercept = 0)+
theme_bw()
grid.arrange(plot1,plot2,plot3,plot4,plot5)
Se puede apreciar cierta linealidad en todas las variables, sin embargo de antemano ya sabemos que el valor que realmente afecta al valor de producción de la papa es la superficie cosechada.
¿Los residuos son normales?
Gráfica de cuantiles
qqnorm(modelo$residuals)
qqline(modelo$residuals)
En la gráfica se puede ver que si son normales, ahora haremos una prueba de normalidad de shapiro-wilk, para comprobar que si sea asi.
Prueba de normalidad Shapiro-Wilk
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.86492, p-value = 0.003427Tanto la gráfica como este test nos muestra que los residuos son normales.
Ahora se probará la variabilidad constante de los residuos (Homocedasticidad).
ggplot(papas, aes(modelo$fitted.values, modelo$residuals))+
geom_point()+
geom_smooth(color="darkred",se=FALSE)+
geom_hline(yintercept = 0)+
theme_bw()
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 5.1179, df = 4, p-value = 0.2754corrplot(cor(dplyr::select(papas, Superficie_Sembrada, Superficie_Cosechada,Superficie_Siniestrada)),
method = "number", tl.col = "black")
Análisis de varianza
vif(lm(Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada, data = papas ))## Superficie_Sembrada Superficie_Cosechada
## 8.157059 8.157059dwt(modelo, alternative="two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.1487266 1.681133 0.194
## Alternative hypothesis: rho != 0Identificacion de posibles valores atípicos o influyentes
papas$studentized_residual <- rstudent(modelo)
ggplot(papas, aes(x=predict(modelo),y=abs(studentized_residual)))+
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución de los residuos studentized",
x = "predicción modelo") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))
which(abs(papas$studentized_residual) > 3)## [1] 20No se observa ningua observacion atipica
summary(influence.measures(modelo))## Potentially influential observations of
## lm(formula = Valor_de_la_Produccion ~ Superficie_Sembrada + Superficie_Cosechada + Superficie_Siniestrada + lluvias + temperaturas, data = papas) :
##
## dfb.1_ dfb.Sp_S dfb.Sp_C dfb.llvs dfb.tmpr dffit cov.r cook.d hat
## 15 -0.05 -0.30 0.27 -0.18 0.07 -0.39 2.43_* 0.03 0.49
## 20 -1.46_* -0.69 1.07_* 0.90 1.34_* 2.36_* 0.01_* 0.45 0.15
## 22 -0.05 0.49 -0.46 -0.17 0.07 0.54 4.49_* 0.06 0.72_*
## 25 1.61_* 0.17 0.00 -1.03_* -1.53_* 2.01_* 3.37_* 0.79 0.73_*influencePlot(modelo)
## StudRes Hat CookD
## 16 -1.8271765 0.1518627 0.10704110
## 20 5.5391152 0.1535683 0.44818077
## 22 0.3343268 0.7204684 0.06029557
## 25 1.2060231 0.7344504 0.78668238Conclusión:
Podemos remarcar que los datos siniestrados se comportan de una manera muy variada y esto se debe a el cambio climático que llega a suceder en ciertas estaciones del año, esto causa confusiones y problemas a los agricultores ya que muchas veces es de imprevisto que una temepratura extrema o no planeada ocurra y eso conduce a la superficie sinistreada de la que es incierto saber su comportamiento,esto con lleva a que no haya una correlación muy natural de esa parte de los datos, por otro lado podemos ver que la superficie cosechada y la sembrada van de la mano ya que como normalmente no se tiene mucha pérdida en los sembradíos estos estarán actuando de la misma manera a lo largo de su periodo, por último cabe mencionar que el valor de producción si es bastante variado ya que pueden ocurrir muchos erxcenarios, unos de ellos es que normalmente se siembran distintas hortalizas en el mismo campo para mantener los nutrientes al nivel y poder aprovechar la tierra mejor.
Bibliografía:
Editor. (2015, 12 noviembre). El poder de. . . La papa. El Poder del Consumidor. https://elpoderdelconsumidor.org/2015/11/el-poder-de-la-papa/
Panorama. (2018, 16 abril). Papa. Panorama AGROPECUARIO. https://panorama-agro.com/?page_id=2547