En este documento reportaré el procedimiento seguido para crear un algoritmo que prediga la fecha de la cosecha tomando en cuenta la variedad, carga y fecha de aplicación de cianamida, tomando en cuenta datos de clima histórico.
Para poder cumplir con el objetivo de predecir la fecha de cosecha en base a variedad, carga, clima y fecha de aplicación de cianamida tenemos que hacer lo siguiente.
Calcular la cantidad de horas/calor sobre 16 grados (solo los positivos) para la data histórica. Sumar el total por el dÃa. Promediar las fechas iguales (ignorando el año).
Predecir la cantidad de horas/sol necesarias (del año especÃfico que tenemos la data, no del histórico), la carga y la variedad.
En base a ese modelo, podemos predecir la cantidad de horas que son necesarias para la madurez dependiendo de la carga y variedad
Buscar en la data histórica cuantos dÃas son necesarios desde la aplicación de la cianamida para lograr esa cantidad predicha de horas sol
Sumar esa cantidad de dÃas a la aplicación de la cianamida para tener la fecha estimada de cosecha
Comenzamos por cargar los datos de Clima (Clima) y de maduración de uvas (Uvas). En ambos casos estamos seleccionando únicamente las variables que vamos a emplear. En cuanto a la temperatura, estamos usando la variable Temperatura Exterior como indicador, y dejando de lado la máxima y la mÃnima que también están disponibles en la base de datos.
La calidad de un modelo predictivo, dependerá de que los insumos de los que se alimenta sean precisos. Vamos a revisar que no hayan cosas extrañas en los datos de la fuente para asegurarnos de la calidad de los datos.
Estos son los datos que tenemos de temperatura cada intervalos de 15 min. Las temperaturas oscuilan entre los 10 y 35 grados, y es notoria la ausencia de registro durante varios periodos antes del 2017.
En la mayor parte de los datos de clima, están dados por intervalos de 15 minutos. Sin embargo, no hay 94 registros en todos los dÃas. En el gráfico vemos en verde los dias con 96 registros y en rojo los dÃas con más o menos registros. Esto es porque hay dÃas con registros faltantes, o dÃas donde se registra la temperatura cada 5 o 10 minutos. Para uniformizar el registro de clima, promediaremos la temperatura a nivel de hora.
Podemos ver que aún hay 15 dÃas con menos registros. Esto es un problema porque esto puede generar vacÃos en el registro del clima promedio, por lo que será necesario imputar estos valores con valores probables.
Después de un proceso de limpieza de datos, ahora tenemos unicamente registros de temperatura promedio a nivel de hora. La falta de estandarización a nivel de intervalo de datos nos impide recuperar la mayor granularidad en los datos. Los registros ausentes fueron completados con el algoritmo mice, dado que no es posible realizar una suma del total del dÃa si el dÃa no cuenta con un registro completo de 24 horas.
El primer paso es definir que indice de calor es el más apropiado. El IC se define como la sumatoria de grados por encima de cierto nivel (usualmente se usa 10C) dentro de cada hora, despues de la aplicación de la cianamida. Esto consistirÃa (aproximadamente en el area pintada de azul en el siguiente gráfico), con la salvedad de que se calcula desde la aplicación de la cianamida hasta la fecha de cosecha ideal.
El IC más apropiado es aquel que tenga la mayor correlacion con la fecha ideal de cosecha.
Debido a como se relacionan los indices de calor con la fecha de cosecha, usaremos el indice de calor en base a 18 grados.
En toda la data historica calcularemos el indice de calor por dÃa. Para esto: 1. Restamos 18 a la temperatura registrada cada 15 minutos 2. Tomamos solo los valores positivos 3. Los sumamos 4. Los promediamos de acuerdo a la fecha, sin tomar en cuenta el año.
El gráfico muestra el nivel promedio en la data histórica de horas calor acumuladas en cada dÃa del primer mes del año.
Ahora, dependiendo de la variedad y la carga, construiremos un modelo que prediga cuanto IC requiere la planta para alcanzar la madurez.
ModHoras = lm(B18 ~ Variedad + Carga,data = Uvas)
ModHoras %>% summary()
##
## Call:
## lm(formula = B18 ~ Variedad + Carga, data = Uvas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2129.5 -796.5 -88.0 842.8 2733.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11576.027 1610.841 7.186 1.03e-08 ***
## VariedadCandy Snaps -4513.927 1828.751 -2.468 0.0179 *
## VariedadCotton Candy -1653.285 1408.325 -1.174 0.2474
## VariedadJack's Salute 494.486 1341.733 0.369 0.7144
## VariedadSweet Celebration -73.319 1399.333 -0.052 0.9585
## VariedadSweet Globe -458.214 1283.715 -0.357 0.7230
## VariedadSweet Jubile -603.359 1794.777 -0.336 0.7385
## VariedadSweet Sapphire 286.842 1406.685 0.204 0.8395
## Carga -1.605 0.360 -4.458 6.54e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1257 on 40 degrees of freedom
## Multiple R-squared: 0.3943, Adjusted R-squared: 0.2732
## F-statistic: 3.255 on 8 and 40 DF, p-value: 0.006
Podemos ver que la cantidad de IC necesaria no es dependiente de manera muy cercana con la variedad y la carga (Solo se explica el 39.43% de la variabilidad en el IC). Además, podemos ver que la carga se relaciona negativamente con la cantidad de horas necesarias, Por cada caja de carga disminuye el IC por 1.61. Esto no corresponde con la información agricultural previa. Esto podrÃa ser un indicador de error en la calidad de los datos. La ecuación completa es: \[B18 = 11576.0 - 4513.9 * VariedadCandy Snaps - 1653.3 * VariedadCotton Candy + 494.5 * VariedadJack's Salute - \\ 73.3 * VariedadSweet Celebration - 458.2 * VariedadSweet Globe - 603.4 * VariedadSweet Jubile + 286.8 * VariedadSweet Sapphire\\ - 1.6 * Carga\]
Una vez que tenemos el modelo de regresión, podemos usar la formula para predecir cuantas horas necesita una planta dependiendo de su variedad y su carga. En el ejemplo podemos ver que podemos esperar que una planta de Candy Snaps con una carga de 3500 cajas necesitará de 1445 de IC para alcanzar el punto optimo de madurez.
predict(ModHoras,newdata = data.frame(Variedad = "Candy Snaps",Carga = 3500))
## 1
## 1445.397
Por ahora, podemos calcular el IC necesario para que la planta madure, pero necesitamos integrar esta información con nuestro registro histórico de clima, de manera que podamos crear una función que reciba como input la fecha de aplicación de la cianamida, calcule cual es el IC que la planta va acumulando, y devuelva la fecha de cosecha estimada.
FechaCosecha("Candy Snaps",3500,"2019-06-23")
## [1] "2019-09-07"
Vamos a calcular que fecha de cosecha hubiera predicho el modelo para los datos que tenemos y luego compararlos con la fecha ideal de cosecha reportada en la base de datos. Depués podemos calcular la diferencia en dÃas entre la fecha real y la fecha estimada en base al algoritmo.
En la figura podemos ver que el algoritmo sistemáticamente subestima la fecha de cosecha (La fecha real siempre es mayor a la estimada).Relación entre fecha real y fecha predicha de cosecha
Si calculamos el tamaños de la desviación, podemos ver que el algorimo falla en promedio por 35.2449 dÃas. El histograma muestra la distribución de los errores del promedio (A), y luego la desviacion de los errores si sumamos los 35 dÃas a la predicción para al menos centrar los errores en 0 (B). En ambos casos los resultados son desalentadores.
Un posible problema es la diferencia entre las horas calor reales acumuladas en el 2019 y las horas calor promedio/históricas
El gráfico explica el problema. La data usada para el modelo es la data del año 2019, pero la data histórica es significativamente diferente de esta. Si fueran iguales, todos los puntos deberÃan caer sobre la lÃnea. Sin embargo los datos reales demuestran que la temperatura históricos es significativamente más alta que la del 2019, y el ángulo de la linea, ni siquiera es paralelo a la de los puntos lo que indica que el problema tampoco se puede solucionar restando una constante a los datos históricos.
Vamos a ver si la diferencia se encuentra en que las temperaturas son demasiado diferentes, o es una cuestión del calculo del IC. Curiosamente, el gráfico muestra las desviaciones en temperatura por hora entre la primera fecha de aplicación de cianamida hasta la última fecha de cosecha. Podemos ver que la diferencia de temperaturas entre el histórico y el real no muestra un sesgo claro. Los puntos rojos muestran las horas en las que la temperatura estuvo por encima en la data histórica, mientras que los rojos muestran las horas en las que la temperatura mfue mayor en el año 2019.