¿Por qué las abejas importan tanto para la humanidad?

Las abejas no solo son extremadamente importantes para los humanos, sino también para el funcionamiento de ecosistemas enteros. Como sabemos, las abejas permiten que las plantas se reproduzcan mediante la polinización. Estas plantas contribuyen al sistema alimentario al alimentar a los animales, además de los humanos, como las aves y los insectos. Si la fuente de alimento para estos animales se redujera o se perdiera por completo, causaría sufrimiento a toda la cadena alimentaria.

Un nuevo estudio publicado en Nature Communications encontró que solo el dos por ciento de las especies de abejas silvestres contribuye con el 80 por ciento de las visitas de polinización de cultivos observadas a nivel mundial. Esto significa que si este pequeño porcentaje de abejas desaparece entonces el 80 por ciento de nuestro sistema agrícola colapsará. Puede parecer increíble, pero sin las abejas, podemos despedirnos de alimentos como manzanas, almendras, naranjas y aguacates. Si bien este es posiblemente el impacto más apremiante de la pérdida de abejas, también hay consecuencias económicas muy reales ya que las abejas contribuyen con alrededor de $ 4.2 billones a la economía global.

¿Por qué las abejas están desapareciendo tan rápidamente?

Existe una variedad de amenazas que enfrenta la población de abejas, incluida la pérdida de hábitat y el cambio climático , pero la amenaza más urgente para las abejas son los pesticidas.

La agricultura ha crecido hasta depender en gran medida de los insecticidas para proteger los cultivos de la destrucción, pero si bien las abejas no amenazan los cultivos, la pulverización de productos químicos a gran escala puede tener un gran impacto en otros insectos, incluidas las moscas y, por supuesto, las abejas. Investigaciones recientes revelaron que muchos pesticidas populares que todavía se usan en los Estados Unidos en realidad alteran la mente de las moscas y las abejas, destruyen su memoria y desechan su ciclo de sueño / vigilia.

Factores que influyen en la muertes de las abejas

En las últimas décadas se estima que hubo un declive de alrededor del 25% de las especies de abejas durante los años de 2006 a 2015 con respecto a los números registrados en el año de 1990, aunque estas interpretaciones deberían ser tomadas con cuidado debido a la naturaleza de los datos inesperados del medio ambiente, se debe mantener una constante revisión sobre ellos puesto que la realidad podría ser peor. Existen muchos factores que han afectado la vida de las abejas e insectos en general, situaciones como la agricultura, cambio climático, deforestación, urbanización y muchos otros pueden ser el motivo por el cual el número de las abejas se encuentra en declive pero de igual manera el aumento de muchas de esas actividades o fenómenos sigue poniendo en riesgo su capacidad de sobrevivencia conforme pasa el tiempo y con ello a presentar las consecuencias de un ecosistema escaso de abejas.

Las abejas no solo son extremadamente importantes para los humanos, sino también para el funcionamiento de ecosistemas enteros. Como sabemos, las abejas permiten que las plantas se reproduzcan mediante la polinización, estas plantas contribuyen al sistema alimentario al alimentar a los animales, además de los humanos, como las aves y los insectos. Si la fuente de alimento para estos animales se redujera o se perdiera por completo, causaría estragos en toda la cadena alimentaria.

¿Cómo es el panorama de las colmenas en México?

Actualmente el número de colmenas en México esta reduciendose con el paso de los años lo que por consecuencia afecta directamente al número de abejas que se tienen en el país y debido a este evento se espera que las poblaciones de abejas tanto en México como en el mundo se encuentren en una situación de peligro, motivo que tiene en alerta al mundo puesto podria ser el inicio de una cadena de desastres naturales de los que no se tendría control.

setwd("~/estadistica aplicada")
library(pacman)
library(readxl)
library(plotrix)
library(ggplot2)
library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:plotrix':
## 
##     rescale

library(RColorBrewer)
p_load("base64enc", "htmltools", "mime", "xfun", "prettydoc", "readr", "ggplot2", "tidyr", "plotly", "readr")

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Bs <- read.csv("BeehivesMX.csv")
Bs %>%
  tail(58) %>%
  ggplot( aes(x=yearc, y=Colmenas)) +
    geom_line(color = "lightblue") +
    geom_point(shape=21, color="orange", fill="yellow", size=3)+
  ggtitle("Evolución de Colmenas en México de 1961 a 2018")

Figura 1. Grafico de número de colmenas en México (1961-2018)

En este gráfico podemos notar como existe un descenso de colmenas a partir del 84 el cual no se ha podido alcanzar desde ese entonces además de mencionar que hubo una epoca despues de 1960 donde hubo un descenso extremo en la cantidad de colmenas, donde uno de los principales motivos a los que se le puede atribuir dicho suceso es a la llamada Revolución Verde, un evento historico en el que la producción agrícola tuvo un avance inmenso con respecto a las anteriores metodologías pero que con ella llegaron los plaguicidas, los cuales afectaron fuertemente a los seres vivos debido al uso desmedido que durante un tiempo se le dio.

A continuacion se vera la evolucion de la exportación de miel que se ha dado anualmente en México, uno de los grandes motivos por los que se realiza la actividad de la apicultura en México

EM <- read.csv("exportamiel.csv")
EM %>%
  tail(26) %>%
  ggplot( aes(x=YEAR, y=MielProd)) +
    geom_line(color = "lightblue") +
    geom_point(shape=21, color="orange", fill="yellow", size=3)+
  ggtitle("Exportaciones de miel en México")

Figura 2. Gráfico de exportacion de miel en México

Con este gráfico podemos ver un panorama más positivo donde México en los años más recientes ha exportado mayores cantidades de miel lo que podria darnos una idea del interes por esta practica pero que a pesar de ello, en parte de la produccion de este tipo solo se hace como un medio de obtencion de bienes y desgraciadamente no estan tan enfocados en un reestablecimiento ecologico a excepcion de algunos.

• ¿CÓMO SE RELACIONA CON LA AGRICULTURA?

En cualquier actividad que se beneficie de las plantas tendra una fuerte relacion con las abejas y su trabajo como polinizadoras debido a que ayudan en gran parte a que los productos que se obtengan sean de buena calidad o en gran cantidad.

Actualmente se utilizan grandes cantidades de colmenas para la polinizacion de huertos dentro de la agricultura, esta se beneficia de las abejas para permitir de esta manera que sus plantas tengan una mejor produccion

Poli <- read.csv("poliniza.csv")
x <- ggplot(data=Poli, aes(x=Delegacion, y=Colmenas)) + 
    geom_bar(stat="identity", position="stack", fill = "orange")
x + coord_flip()

Figura 3. Gráfico de colmenas en usos agrícolas en 2008

Con esto podemos ver como existe una gran cantidad de colmenas en esta industria, lo que nos permite reconocer que incluso con el gran trabajo que realizan las abejas por el medio ambiente hay maneras en las que ambas industrias pueden coexistir y trabajar juntas para la obtención de un bien comun algo de gran utilidad.

¿Qué sucede?

En los recientes años se ha visto un gran descenso en la cantidad de abejas y respectivamente de colmenas alrededor de todo el mundo y México, considerandose una alerta para una posible extinción de estas mismas, lo que nos llevaria a un desequilibrio natural de gran magnitud del cual dificilmente podriamos recuperarnos, si bien puede ser producto de multiples factores podemos decir que la producción de recursos economicos es uno de los grandes motivos por el cual la naturaleza de nuestro planeta se ve afectada, actividades primarias como la agrgicultura son una de ellas puesto que sin importar el daño que se le realice a la vida silvestre donde se cultive deciden no ver por ellas y seguir con sus actividades.

Los plaguicidas son uno de los motivos más grandes por los que insectos como las abejas mueren en regiones del país con mucha agricultura, el solo hecho de llevar una produccion de plaguicidas es suficiente para generar un impacto ambiental, desde su consumo hasta su desperdicio, por lo que hemos realizado un estudio de la producción de plaguicidas en toneladas con respecto a las colmenas desde 1994 hasta 2017.

Datos para análisis de correlación lineal

setwd("~/estadistica aplicada")
library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:plotly':
## 
##     select

library(ggplot2)
library(pacman)
library(readxl)
PB <- read_excel("plaguibes.xlsx")

xfun::embed_file("plaguibes.xlsx")

Download plaguibes.xlsx

Diagrama de dispersión

ggplot(data = PB, aes(x = tonPlagui, y = Colmenas)) + 
  geom_point(colour = "red4") +
  ggtitle("Diagrama de dispersión") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

El diagrama de dispersion nos permite ver la manera en la que se distribuyen puntualmente cada dato lo que nos permite ver si existe alguna aparente relacion entre ellas.

Análisis de normalidad

par(mfrow = c(1, 2))
hist(PB$tonPlagui, breaks = 10, main = "", xlab = "Plaguicidas", border = "darkred")
hist(PB$Colmenas, breaks = 10, main = "", xlab = "Colmenas",
     border = "blue")

qqnorm(PB$tonPlagui, main = "Plaguicidas", col = "darkred")
qqline(PB$tonPlagui)

qqnorm(PB$Colmenas, main = "Colmenas", col = "blue")
qqline(PB$Colmenas)

par(mfrow = c(1,1))

Significancia

Se prosigue a realizar un test de hipotesis de análisis de normalidad.

shapiro.test(PB$tonPlagui)

## 
##  Shapiro-Wilk normality test
## 
## data:  PB$tonPlagui
## W = 0.91665, p-value = 0.04928

shapiro.test(PB$Colmenas)

## 
##  Shapiro-Wilk normality test
## 
## data:  PB$Colmenas
## W = 0.86951, p-value = 0.005147

• Representación Grafica

par(mfrow = c(1, 2))
hist(log10(PB$Colmenas), breaks = 10, main = "", xlab = "Log10(Colmenas)",
     border = "blue")
qqnorm(log10(PB$Colmenas), main = "", col = "blue")
qqline(log10(PB$Colmenas))

par(mfrow = c(1, 1))
shapiro.test(log10(PB$Colmenas))

## 
##  Shapiro-Wilk normality test
## 
## data:  log10(PB$Colmenas)
## W = 0.89623, p-value = 0.01792

Homoestaticidad

ggplot(data = PB, aes(x = tonPlagui, y = Colmenas)) + 
  geom_point(colour = "red4") +
  geom_segment(aes(x = 10000, y = 1790000, xend = 11000, yend = 2300000),linetype="dashed") +
  geom_segment(aes(x = 10000, y = 1700000, xend = 38000, yend = 1720000),linetype="dashed") +
  ggtitle("Diagrama de dispersión") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

Cálculo de correlación

cor(x = PB$tonPlagui, y = log10(PB$Colmenas), method = "pearson")

## [1] -0.3362168

cor(x = PB$tonPlagui, y = log10(PB$tonPlagui), method = "spearman")

## [1] 1

Significancia de correlación

cor.test(x = PB$tonPlagui,
         y = log10(PB$Colmenas), 
         alternative = "two.sided",
         conf.level  = 0.95,
         method      = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  PB$tonPlagui and log10(PB$Colmenas)
## t = -1.6745, df = 22, p-value = 0.1082
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.65128105  0.07772117
## sample estimates:
##        cor 
## -0.3362168

cor.test(x = PB$tonPlagui,
         y = log10(PB$Colmenas),
         alternative = "two.sided",
         conf.level  = 0.95,
         method      = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  PB$tonPlagui and log10(PB$Colmenas)
## S = 2956, p-value = 0.1762
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2852174

Coeficiente de determinació R^2

R2_pearson <- cor(x = PB$tonPlagui,
                  y = log10(PB$Colmenas),
                  method = "pearson")
R2_pearson <- R2_pearson^2
R2_pearson

## [1] 0.1130417

R2_spearman <- cor(x = PB$tonPlagui,
                  y = log10(PB$Colmenas),
                   method = "spearman")
R2_spearman <- R2_spearman^2
R2_spearman

## [1] 0.08134896

Con los datos obtenidos podemos ver que puede no existir un relación directa entre el numero de colmenas con respecto a la producción de plaguicidas, donde a base de análisis podriamos decir que puede ser debido a los múltiples factores que también aportan dentro del descenso de las abejas e incluso con el hecho de ser solo una base de producción y no de consumo como tal pero que estando en México podria ser conveniente realizar una muestra con solo las regiones agrícolas del país.

Datos para regresión lineal múltiple

setwd("~/estadistica aplicada")
library(pacman)
library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:scales':
## 
##     alpha, rescale

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

## The following object is masked from 'package:plotrix':
## 
##     rescale

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()
## )

Análisis para regresión lineal múltiple

• 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.000

Una 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   ggplot2

ggpairs(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`.

Generación 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.00896

Selección de 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-02

El 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.00896

Es 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-02

Cada 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)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

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.127

Se 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.numeric

bptest(modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 11.096, df = 4, p-value = 0.02551

No hay evidencias de homocedasticidad

Matríz de correlación entre predictores

library(corrplot)

## corrplot 0.84 loaded

corrplot(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':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

vif(modelo)

##       Deforestacion           TonPlagui           Tempmedia Densidadpoblacional 
##            1.325878            2.920699            2.095663            4.976713

Los 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.188
##  Alternative hypothesis: rho != 0

No hay evidencia de autocorrelación

Identificació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 
## 17

Se 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.39

En 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.6632014