APICULTURA - Regresion Multiple
INTRODUCCIÓN
La apicultura es uno de los sectores más importantes dentro de las actividades economicas puesto que es de la unica manera que se puede obtener miel y cera de abeja, además de que se basa en la preservación de la población de abejas y que incluso tiene un impacto positivo dentro del medio ambiente, beneficiando de gran manera las conidiciones en la que se encuentra la tierra.
Debido al consumo de los medios naturales, cambio climatico y emisiones de plaguicidas se esta en decadencia la población de abejas actualmente, llegando a números criticos que amenazan a toda la diversidad ecologica por lo que es necesario tomarle importancia a esta situación.
OBJETIVO
Realizar un analisis de Regresión Lineal Multiple con respecto a la situacion de colmenas de abejas en México con algunos factores que afecta a su preservación.
REGRESION LINEAL MÚLTIPLE
La regresión lineal múltiple permite generar un modelo lineal en el que el valor de la variable dependiente o respuesta (Y ) se determina a partir de un conjunto de variables independientes llamadas predictores (X1 , X2 , X3 …). Es una extensión de la regresión lineal simple, por lo que es fundamental comprender esta última. Los modelos de regresión múltiple pueden emplearse para predecir el valor de la variable dependiente o para evaluar la influencia que tienen los predictores sobre ella (esto último se debe que analizar con cautela para no malinterpretar causa-efecto).
Los modelos lineales múltiples siguen la siguiente ecuación:
\(Y_{i}=(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i}+\cdots+\beta_{n}X_{ni})+e_{i}\)
β0 es la ordenada en el origen, el valor de la variable dependiente Y cuando todos los predictores son cero.
βi : es el efecto promedio que tiene el incremento en una unidad de la variable predictora Xi sobre la variable dependiente Y , manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.
ei: es el residuo o error, la diferencia entre el valor observado y el estimado por el modelo.
Es importante tener en cuenta que la magnitud de cada coeficiente parcial de regresión depende de las unidades en las que se mida la variable predictora a la que corresponde, por lo que su magnitud no está asociada con la importancia de cada predictor. Para poder determinar qué impacto tienen en el modelo cada una de las variables, se emplean los coeficientes parciales estandarizados, que se obtienen al estandarizar (sustraer la media y dividir entre la desviación estándar) las variables predictoras previo ajusteS
Datos
setwd("~/estadistica aplicada")
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()
## )Analisis
Relacion entre las variables
Esta información es crítica a la hora de identificar cuáles pueden ser los mejores predictores para el modelo, qué variables presentan relaciones de tipo no lineal (por lo que no pueden ser incluidas) y para identificar colinialidad entre predictores
- Matriz de Coeficientes de correlacion
round( cor( x = datos, method = "pearson"), 3)## TonPlagui Colmenas Deforestacion Tempmedia
## TonPlagui 1.000 0.420 0.370 0.284
## Colmenas 0.420 1.000 0.321 0.422
## Deforestacion 0.370 0.321 1.000 0.318
## Tempmedia 0.284 0.422 0.318 1.000
## Densidadpoblacional 0.763 0.735 0.496 0.647
## Densidadpoblacional
## TonPlagui 0.763
## Colmenas 0.735
## Deforestacion 0.496
## Tempmedia 0.647
## Densidadpoblacional 1.000Una vez calculado el coeficiente de correlación de cada una de las variables, podemos concluir lo siguiente para las variables más importantes:
Podemos obsrvar que la variable colmenas tiene mayor correlacion con la densidad poblacional. Más sin embargo, esta variable tiene una baja correlación con la deforestacion,los plaguisidas y tempertura.
- Distribucion de los datos
library(psych)
multi.hist(x = datos, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
main = "")
- GGally
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(datos, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Generacion del modelo
Ahora una vez que entendemos la forma en la cual se relacionan las variables, podemos empezar a experimentar con la generacion del modelo
modelo <- lm(Colmenas ~ Deforestacion + TonPlagui + Tempmedia + Densidadpoblacional , data = datos )
summary(modelo)##
## Call:
## lm(formula = Colmenas ~ Deforestacion + TonPlagui + Tempmedia +
## Densidadpoblacional, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -140422 -29068 -11899 21169 147205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.932e+04 4.432e+05 0.179 0.86073
## Deforestacion -8.448e-02 2.604e-01 -0.324 0.75078
## TonPlagui -3.950e+00 2.451e+00 -1.611 0.13108
## Tempmedia -2.774e+05 2.629e+05 -1.055 0.31065
## Densidadpoblacional 1.863e-02 5.479e-03 3.401 0.00474 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82410 on 13 degrees of freedom
## Multiple R-squared: 0.6226, Adjusted R-squared: 0.5065
## F-statistic: 5.361 on 4 and 13 DF, p-value: 0.00896Selección de los mejores predictores
step(object = modelo, direction = "both", trace = 1)## Start: AIC=411.64
## Colmenas ~ Deforestacion + TonPlagui + Tempmedia + Densidadpoblacional
##
## Df Sum of Sq RSS AIC
## - Deforestacion 1 7.1483e+08 8.9009e+10 409.79
## - Tempmedia 1 7.5590e+09 9.5853e+10 411.12
## <none> 8.8294e+10 411.64
## - TonPlagui 1 1.7637e+10 1.0593e+11 412.92
## - Densidadpoblacional 1 7.8544e+10 1.6684e+11 421.10
##
## Step: AIC=409.79
## Colmenas ~ TonPlagui + Tempmedia + Densidadpoblacional
##
## Df Sum of Sq RSS AIC
## - Tempmedia 1 7.5116e+09 9.6520e+10 409.25
## <none> 8.9009e+10 409.79
## - TonPlagui 1 1.7514e+10 1.0652e+11 411.02
## + Deforestacion 1 7.1483e+08 8.8294e+10 411.64
## - Densidadpoblacional 1 8.0436e+10 1.6944e+11 419.38
##
## Step: AIC=409.25
## Colmenas ~ TonPlagui + Densidadpoblacional
##
## Df Sum of Sq RSS AIC
## - TonPlagui 1 1.1140e+10 1.0766e+11 409.21
## <none> 9.6520e+10 409.25
## + Tempmedia 1 7.5116e+09 8.9009e+10 409.79
## + Deforestacion 1 6.6742e+08 9.5853e+10 411.12
## - Densidadpoblacional 1 9.6219e+10 1.9274e+11 419.70
##
## Step: AIC=409.21
## Colmenas ~ Densidadpoblacional
##
## Df Sum of Sq RSS AIC
## <none> 1.0766e+11 409.21
## + TonPlagui 1 1.1140e+10 9.6520e+10 409.25
## + Tempmedia 1 1.1377e+09 1.0652e+11 411.02
## + Deforestacion 1 5.8675e+08 1.0707e+11 411.12
## - Densidadpoblacional 1 1.2629e+11 2.3395e+11 421.18##
## Call:
## lm(formula = Colmenas ~ Densidadpoblacional, data = datos)
##
## Coefficients:
## (Intercept) Densidadpoblacional
## 6.559e+05 1.059e-02El mejor modelo resultante del proceso de selección ha sido:
modelo <- (lm(formula = Colmenas ~ Deforestacion + TonPlagui + Tempmedia + Densidadpoblacional, data = datos))
summary(modelo)##
## Call:
## lm(formula = Colmenas ~ Deforestacion + TonPlagui + Tempmedia +
## Densidadpoblacional, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -140422 -29068 -11899 21169 147205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.932e+04 4.432e+05 0.179 0.86073
## Deforestacion -8.448e-02 2.604e-01 -0.324 0.75078
## TonPlagui -3.950e+00 2.451e+00 -1.611 0.13108
## Tempmedia -2.774e+05 2.629e+05 -1.055 0.31065
## Densidadpoblacional 1.863e-02 5.479e-03 3.401 0.00474 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82410 on 13 degrees of freedom
## Multiple R-squared: 0.6226, Adjusted R-squared: 0.5065
## F-statistic: 5.361 on 4 and 13 DF, p-value: 0.00896Es recomendable mostrar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:
confint(lm(formula = Colmenas ~ Deforestacion + TonPlagui + Tempmedia +
Densidadpoblacional, data = datos))## 2.5 % 97.5 %
## (Intercept) -8.781692e+05 1.036803e+06
## Deforestacion -6.470481e-01 4.780882e-01
## TonPlagui -9.244529e+00 1.345340e+00
## Tempmedia -8.454025e+05 2.906406e+05
## Densidadpoblacional 6.795538e-03 3.046877e-02Cada una de las pendientes de un modelo de regresión lineal múltiple (coeficientes parciales de regresión de los predictores) se define del siguiente modo: Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable (Y) varía en promedio tantas unidades como indica la pendiente.
Validación de condiciones para la regresión múltiple lineal
Relación lineal entre los predictores numéricos y la variable respuesta:
Esta condición se puede validar bien mediante diagramas de dispersión entre la variable dependiente y cada uno de los predictores (como se ha hecho en el análisis preliminar) o con diagramas de dispersión entre cada uno de los predictores y los residuos del modelo. Si la relación es lineal, los residuos deben de distribuirse aleatoriamente en torno a 0 con una variabilidad constante a lo largo del eje X. Esta última opción suele ser más indicada ya que permite identificar posibles datos atípicos.
library(ggplot2)
library(gridExtra)
datos <- read_csv("abejas.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## TonPlagui = col_double(),
## Colmenas = col_double(),
## Deforestacion = col_double(),
## Tempmedia = col_double(),
## Densidadpoblacional = col_double()
## )plot1 <- ggplot(data = datos, aes(TonPlagui, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(Tempmedia, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(Densidadpoblacional, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(Deforestacion, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot5 <- ggplot(data = datos, aes(Colmenas, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4, plot5)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Distribucion normal de residuos
qqnorm(modelo$residuals)
qqline(modelo$residuals)
Shapiro Test
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.91957, p-value = 0.127Se confirma la normalidad.
Variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo, los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 11.096, df = 4, p-value = 0.02551No hay evidencias de homocedasticidad
Matriz de correlación entre predictores
library(corrplot)## corrplot 0.84 loadedcorrplot(cor(dplyr::select(datos, TonPlagui, Colmenas, Deforestacion, Tempmedia, Densidadpoblacional)),
method = "number", tl.col = "black")
Análisis de inflación de varianza (VIF)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:psych':
##
## logitvif(modelo)## Deforestacion TonPlagui Tempmedia Densidadpoblacional
## 1.325878 2.920699 2.095663 4.976713Los predictores muestran una correlación lineal e inflación de varianza algo alta.
Autocorrelación
library(car)
dwt(modelo, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 -0.1259029 1.786772 0.184
## Alternative hypothesis: rho != 0No hay evidencia de autocorrelación
IDENTICICACIÓN DE POSIBLES VALORES ATÍPICOS O INFLUYENTES
datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
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(datos$studentized_residual) > 3)## 17
## 17Se indentifico un valor atípico
summary(influence.measures(modelo))## Potentially influential observations of
## lm(formula = Colmenas ~ Deforestacion + TonPlagui + Tempmedia + Densidadpoblacional, data = datos) :
##
## dfb.1_ dfb.Dfrs dfb.TnPl dfb.Tmpm dfb.Dnsd dffit cov.r cook.d hat
## 2 0.01 0.03 -0.08 0.00 0.01 0.13 2.25_* 0.00 0.35
## 17 1.75_* 2.61_* 0.51 -0.25 -1.55_* -3.29_* 0.14 1.27_* 0.51
## 18 -1.78_* -1.31_* -0.76 -0.82 1.73_* 2.26_* 0.19 0.66 0.39En la tabla generada se recogen las observaciones que son significativamente influyentes en al menos uno de los predictores (una columna para cada predictor). Las tres últimas columnas son 3 medidas distintas para cuantificar la influencia. A modo de guía se pueden considerar excesivamente influyentes aquellas observaciones para las que:
Leverages (hat): Se consideran observaciones influyentes aquellas cuyos valores hat superen 2.5((p+1)/n), siendo p el número de predictores y n el número de observaciones. Distancia Cook (cook.d): Se consideran influyentes valores superiores a 1. La visualización gráfica de las influencias se obtiene del siguiente modo:
influencePlot(modelo)
## StudRes Hat CookD
## 16 -1.143998 0.4603283 0.2180847
## 17 -3.196245 0.5143505 1.2662638
## 18 2.836610 0.3885622 0.6632014CONCLUSIÓN
El modelo de Regresion lineal múltiple es una herramienta que nos sirve dentro del analisis de dependecia y correlacion de algunas variables con respecto a un tema en especifico puesto que nos brinda la manera de comprender el comportamiento de varios datos dentro de un solo esquema de comparación, como lo fue en el caso de las colmenas con cada una de las variables que se vieron, y que con ello podemos llegar a una conclusión acerca de la relación que pueda llegar a exisitir entre estas.