Resumen

Los datos inflados con cero de los experimentos de campo pueden ser problemáticos, ya que estos datos requieren el uso de modelos estadísticos específicos durante el proceso de análisis. Este estudio utilizó el modelo Binomial Negativo Cero Inflado (ZINB) con las funciones de enlace log y logístico. En este estudio se ajusta el número de huevos en la región de xxxx en el periodo 2018 a 2020 en función de variables climáticas.Se realizó una selección de variable, a partir del Análisis de Componentes Principales (\(ACP\)).

Introducción

El psílido asiático de los cítricos Diaphorina citri Kuwayama es importante por ser el vector de la bacteria Candidatus Liberibacter, agente causal de la Huanglongbing HLB o enverdecimiento de los cítricos “greening”, enfermedad devastadora de la citricultura mundial. El insecto causa daños directos, al alimentarse de la planta e indirectos cuando transmite la bacteria causante de la enfermedad; afecta zonas de crecimiento de la planta ocasionando la caída y deformación de hojas, clorosis, necrosis, enanismo y disminución de la calidad de la fruta; las ninfas y los adultos ocasionan daños al extraer gran cantidad de floema de las hojas y los pecíolos; además, poblaciones altas del insecto causan la caída de flores y frutos pequeños [2]. La bacteria, se caracteriza como proteobacteria del genero Candidatus Liberibacter para la enfermedad Huanglongbing, enfermedad generada por Candidatus Liberibacter asiaticus, Candidatus Liberibacter africanus y Candidatus Liberibacter americus. La enfermedad tiene alta incidencia en Asia y África donde se señala como factor limitante para la producción de cítricos. La bacteria es transmitida por el vector Diaphorina citri . En algunos casos, los cultivos pueden contener grandes poblaciones de D. citri. Por tanto, un conjunto de datos que abarque la variable de incidencia de plantas enfermas puede contener muchos ceros. (Esta parte es del paper )

Material y métodos

Suponga que para cada observación hay dos casos posibles. Si ocurre el caso 1, la cuenta es cero. Sin embargo, si ocurre el caso 2, los recuentos (incluidos los ceros) se generan de acuerdo con el modelo Binomial Negativo. Suponga que el caso 1 ocurre con probabilidad \(\pi\) y el caso 2 ocure con probabilidad \((1 - \pi)\). Por lo tanto, la probabilidad de la variable \(y_{i}\) que se distribuye Binomial Negativa Inflada con Ceros por sus siglas en ingles \(ZINB\) se puede expresar como:

\[\begin{equation} Pr\left(y_{i} = j \right) = \begin{cases} \pi_{i} + (1 - \pi_{i})\cdot g(y_{i} = 0) & si \hspace{0.5cm} j=0 \\ (1 - \pi_{i})\cdot g(y_{i}) & si \hspace{0.5cm} j>0 \end{cases} \end{equation}\]

donde \(\pi_{i}\) es la función de enlace logística y \(g(y_{i})\) es la distribución Binomial Negativa expresada como:

\[\begin{equation} g(y_{i}) = Pr(y = y_{i}/\mu_{i}, \alpha) = \frac{\Gamma (y_{i} + \alpha^{-1})}{\Gamma(\alpha^{-1})\cdot \Gamma(y_{i}+1)}\cdot \left[\frac{1}{1 + \alpha \mu_{i}} \right]^{\alpha^{-1}} \cdot \left[\frac{\alpha \mu_{i}}{1+\alpha \mu_{i}} \right]^{y_{i}} \end{equation}\]

El componente Binomial Negativo puede incluir un tiempo de exposición \(t\) y un conjunto de \(k\) variables regresoras. La expresión que relaciona las cantidades es:

\[\begin{equation} \mu_{i} = exp \left(ln(t_{i}) + \beta_{1}x_{1i} + \beta_{2}x_{2i} + \ldots + \beta_{k}x_{ki} \right) \end{equation}\]

La función de enlace logístico \(\pi_{i}\) esta dado por:

\[\begin{equation} \pi_{i} = \frac{\lambda_{i}}{1 + \lambda_{i} } \end{equation}\]

donde,

\[\begin{equation} \lambda_{i} = exp\left( ln(t_{i}) + \gamma_{1} Z_{1i} + \gamma_{2} Z_{2i} + \ldots + \gamma_{m} Z_{mi} \right) \end{equation}\]

El componente logístico incluye un tiempo \(t\) de exposición y \(m\) variables regresoras.

Los coeficientes de regresión son estimados usando el método de Máxima Verosimilitud (MVM). El logaritmo de la función de verosimilitud es:

\[\mathcal{L} = L_{1} + L_{2} + L_{3} + L_{4}\] donde,

\[L_{1} = \sum_{{i:y_{i}=0 }} ln\left( \lambda_{i} + (1+\alpha \mu_{i})^{-\alpha^{-1}} \right)\] \[L_{2} = \sum_{{i:y_{i}>0 }} \sum_{{j=0 }}^{y_{i}-1} ln\left(j+\alpha^{-1} \right)\] \[L_{3} = \sum_{{i:y_{i}>0 }} \left( -ln(y_{i}!) - (y_{i}+\alpha^{-1})ln(1+\alpha \mu_{i})+y_{i}ln(\alpha) + y_{i}ln(\mu_{i}) \right)\] \[L_{4} = \sum_{i=1}^{n} ln(1+\lambda_{i})\]

Resultados

1. Preparación de los datos

En esta sección se organizan los datos para análisis. Las variables de respuesta como (brotes evaluados, número de huevos, ninfas, adultdos, número de huevos por brotes y número de ninfas por brote) se totalizan por semana y las variables climáticas se promedian por semana.

BD_Campo_16_Citricos_F2_2017_2020 <- openxlsx::read.xlsx("BD_Campo_16_Citricos_F2_2017-2020.xlsx",
                                                         sheet = 1)


Fecha <- paste0(BD_Campo_16_Citricos_F2_2017_2020$DIA, "-",
                BD_Campo_16_Citricos_F2_2017_2020$MES, "-", 
                BD_Campo_16_Citricos_F2_2017_2020$AÑO)

Fecha  <- as.Date(Fecha , "%d-%m-%Y")

BD_Campo_16_Citricos_F2_2017_2020 <- cbind.data.frame(Fecha, BD_Campo_16_Citricos_F2_2017_2020)


BD_Campo_16_Citricos_F2_2017_2020$BROTES.EVALUADO[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$HUEVOS[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$NINFAS[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$ADULTOS[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$`N°.DE.HUEVOS.X.BROTES`[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$`N°.DE.NINFAS.X.BROTES`[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$`N°.DE.ADULTOS.X.BROTES`[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0
BD_Campo_16_Citricos_F2_2017_2020$TDIAPHORINA[BD_Campo_16_Citricos_F2_2017_2020$TOTAL.DE.BROTES == 0] <- 0


BD_Campo_16_Citricos_F2_2017_2020 <- BD_Campo_16_Citricos_F2_2017_2020 %>%
                                     dplyr::select(-Tempint15, -HRint15, -HRint10, -Tempint10)

#Filtrar cultivo de naranja
BD_Naranja <- BD_Campo_16_Citricos_F2_2017_2020 %>% 
              filter(AÑO == 2018) %>% 
              filter(CULTIVAR == "Naranja") %>% 
              mutate(Semana = week(Fecha))



#Calcular el promedio por semana
BD_Naranja_Prom_Var_Climaticas <- BD_Naranja  %>%
                                  dplyr::group_by(Semana) %>%
                                  dplyr::select(-DIA, -MES, -AÑO, -MUESTREO, -CULTIVAR, -VARIEDAD,  -PORTAINJERTO,  -REPETICION,
                                                 -"BROTES.EVALUADOS", -NINFAS, -ADULTOS, -"N°.DE.HUEVOS.X.BROTES", 
                                                -"N°.DE.NINFAS.X.BROTES", -"N°.DE.ADULTOS.X.BROTES", -TDIAPHORINA ) %>%
                                  summarise_all(funs(mean(., na.rm = TRUE)))

BD_Naranja$NINFAS <- as.numeric(BD_Naranja$NINFAS)
BD_Naranja$ADULTOS <- as.numeric(BD_Naranja$ADULTOS)

BD_Naranja_Sum_Var_Respuesta <- BD_Naranja  %>%
  dplyr::group_by(Semana) %>%
  dplyr::select("BROTES.EVALUADOS", NINFAS, ADULTOS, "N°.DE.HUEVOS.X.BROTES", 
                "N°.DE.NINFAS.X.BROTES", "N°.DE.ADULTOS.X.BROTES", TDIAPHORINA ) %>%
  summarise_all(funs(sum(., na.rm = TRUE)))

BD_Naranja_Analisis <- cbind.data.frame(BD_Naranja_Sum_Var_Respuesta,
                                        BD_Naranja_Prom_Var_Climaticas[,-1])

2. Análisis Exploratorio de Datos

hist(BD_Naranja_Analisis$HUEVOS, col = 2, main = "Histrograma del número de huevos por semana",
     xlab = "Número de huevos",
     ylab = "Frecuencias")

ggplot(BD_Naranja_Analisis , aes(x = Fecha, y = HUEVOS)) +
  geom_line() + 
  xlab("Semanas")+
  ylab("Total de huevos")

#Promedio del brotes

ggplot(BD_Naranja_Analisis, aes( Fecha)) + 
  geom_line(aes(y = TOTAL.DE.BROTES, colour = "Total de brotes")) + 
  geom_line(aes(y = HUEVOS, colour = "Total de huevos"))+
  xlab("Semanas")+
  ylab("Cantidad")

2.1. Correlación entre las variables

Covariables <- BD_Naranja_Analisis %>%
               dplyr::select(-Fecha, -Semana, -CODIGO,
                             -TOTAL.DE.BROTES,-BROTES.EVALUADO,
                             -BROTES.EVALUADOS, -HUEVOS, -NINFAS,                
                             -ADULTOS, -"N°.DE.HUEVOS.X.BROTES", -"N°.DE.NINFAS.X.BROTES", 
                             -"N°.DE.ADULTOS.X.BROTES", -TDIAPHORINA)  

Covariables <- na.omit(Covariables)



correlations <- cor(Covariables)
corrplot(correlations, tl.col = "gray", method="circle", order = "hclust")

Entre las variables climáticas se observan altas correlaciones entre las variables medidas en diferentes instantes de tiempos, por ejemplo las variables velocidad del viento mínima ,velocidad de viento máxima, evapotranspiración, radiación solar, energía solar, temperatura, temperatura externa, temperatura máxima, temperatura mínima, humedad relativa mínima, humedad relativa máxima, humedad relativa externa, temperatura mínima y lluvia medidas en 10 y 15 días presentan una alta correlación, este comportamiento entre covariables de un modelo generan multicolinedad.

2.2. Análisis de Componentes Principales sobre las variables climáticas

data_PCA <- princomp(Covariables, cor=T)
score <- data_PCA$scores



# Biplot como herramienta exploratoria:
fviz_pca_biplot(data_PCA,
                col.var = "contrib",
                invisible = "ind",
                habillage ="none",
                geom = "text", labelsize=4) + theme_minimal()

El ACP sobre las variables climáticas mostró agrupaciones entre las variables climáticas medidas en diferentes instantes (10 y 15 días). De acuerdo a esto, se procede a dejar las variables Temperatura, Humedad Relativa, Velocidad del Viento, LLuvia, Radiación Solar y Evapo Transpiración medida a los 15 días.

3. Ajuste del modelo \(ZINB\)

library(pscl)

BD_Naranja_Analisis$HuevosEnteros <- round(BD_Naranja_Analisis$HUEVOS,0)

Modelo_3_IC = zeroinfl(HuevosEnteros ~  LLUVIA15, data = BD_Naranja_Analisis, dist = "negbin")
summary(Modelo_3_IC)
## 
## Call:
## zeroinfl(formula = HuevosEnteros ~ LLUVIA15, data = BD_Naranja_Analisis, 
##     dist = "negbin")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0818 -0.6375 -0.5645  0.1213  3.9187 
## 
## Count model coefficients (negbin with log link):
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 0.660238   0.330218   1.999  0.04556 * 
## LLUVIA15    0.006994   0.002210   3.164  0.00156 **
## Log(theta)  3.016341   6.048093   0.499  0.61797   
## 
## Zero-inflation model coefficients (binomial with logit link):
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.496588   0.634123   0.783    0.434
## LLUVIA15    -0.011645   0.009645  -1.207    0.227
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 20.4165 
## Number of iterations in BFGS optimization: 13 
## Log-likelihood: -36.94 on 5 Df
BD_Naranja_Analisis$Ajuste <- predict(Modelo_3_IC, type="response")
BD_Naranja_Analisis <- BD_Naranja_Analisis[with(BD_Naranja_Analisis, order(Semana)), ]
#Ajuste del modelo
ggplot(BD_Naranja_Analisis, aes(x = Semana, y = Ajuste)) +
  geom_line(aes(y = HuevosEnteros, color = "Datos observados"), size = 1.5 )   +
  geom_line(size = 1, aes(color = "Datos estimados"), size = 1.5 ) +
  labs(x = "", y = "Número total de huevos")
## Warning: Duplicated aesthetics after name standardisation: size