Avocados

Actividad 2.1 Aguacate

Contexto

La base de datos avocado.csv contiene datos sobre la venta de aguacates en Estados Unidos, lo cual incluye datos como precios, regiones, volumen de ventas, entre otros. La información que se busca obtener de esta base de datos es analizar el comportamiento de precios y ventas, para posteriormente obtener un pronóstico de ventas.

Análisis

Lectura de la base de datos

av <- read.csv("C:\\Users\\art191127\\Downloads\\avocado.csv")

Carga de librerías

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1     ✔ purrr   1.0.1
## ✔ tibble  3.1.8     ✔ dplyr   1.1.0
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.4     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(lubridate)

## 
## Attaching package: 'lubridate'
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(scales)

## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

options("install.lock"=FALSE)
library(ggmap)

## ℹ Google's Terms of Service: <]8;;https://mapsplatform.google.comhttps://mapsplatform.google.com]8;;>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

library(dplyr)

avocados <- av
names(avocados)

##  [1] "Date"         "AveragePrice" "Total.Volume" "X4046"        "X4225"       
##  [6] "X4770"        "Total.Bags"   "Small.Bags"   "Large.Bags"   "XLarge.Bags" 
## [11] "type"         "year"         "region"

Resumen de la base de datos aún sin modificar.

names(av)[names(av) == "ï..Date"] <- "Date"
names(avocados)[names(avocados) == "ï..Date"] <- "Date"

avocados$Date <- as.Date(avocados$Date, "%d/%m/%Y")
avocados$region <- as.factor(avocados$region)
avocados$type <- as.factor(avocados$type)
summary(avocados)

##       Date             AveragePrice    Total.Volume          X4046         
##  Min.   :2015-01-04   Min.   :0.440   Min.   :      85   Min.   :       0  
##  1st Qu.:2015-10-25   1st Qu.:1.100   1st Qu.:   10839   1st Qu.:     854  
##  Median :2016-08-14   Median :1.370   Median :  107377   Median :    8645  
##  Mean   :2016-08-13   Mean   :1.406   Mean   :  850644   Mean   :  293008  
##  3rd Qu.:2017-06-04   3rd Qu.:1.660   3rd Qu.:  432962   3rd Qu.:  111020  
##  Max.   :2018-03-25   Max.   :3.250   Max.   :62505647   Max.   :22743616  
##                                                                            
##      X4225              X4770           Total.Bags         Small.Bags      
##  Min.   :       0   Min.   :      0   Min.   :       0   Min.   :       0  
##  1st Qu.:    3009   1st Qu.:      0   1st Qu.:    5089   1st Qu.:    2849  
##  Median :   29061   Median :    185   Median :   39744   Median :   26363  
##  Mean   :  295155   Mean   :  22840   Mean   :  239639   Mean   :  182195  
##  3rd Qu.:  150207   3rd Qu.:   6243   3rd Qu.:  110783   3rd Qu.:   83338  
##  Max.   :20470573   Max.   :2546439   Max.   :19373134   Max.   :13384587  
##                                                                            
##    Large.Bags       XLarge.Bags                 type           year     
##  Min.   :      0   Min.   :     0.0   conventional:9126   Min.   :2015  
##  1st Qu.:    127   1st Qu.:     0.0   organic     :9123   1st Qu.:2015  
##  Median :   2648   Median :     0.0                       Median :2016  
##  Mean   :  54338   Mean   :  3106.4                       Mean   :2016  
##  3rd Qu.:  22029   3rd Qu.:   132.5                       3rd Qu.:2017  
##  Max.   :5719097   Max.   :551693.7                       Max.   :2018  
##                                                                         
##                  region     
##  Albany             :  338  
##  Atlanta            :  338  
##  BaltimoreWashington:  338  
##  Boise              :  338  
##  Boston             :  338  
##  BuffaloRochester   :  338  
##  (Other)            :16221

Lo primero que notamos en esta base de datos con dimensiones de 13 columnas con 18,249 renglones es que tiene 3 columnas que carecen de especificación, las cuales son X4046, X42425 y X4770, por lo que se investigo en Kaggle y se descubrió que se refieren a tipos de aguacates, por lo que en orden se denominan como: 1. Small Hass, 2. Large Hass y 3. Extra Large Hass. Motivo por el cual uno de los primeros pasos será renombrar estas columnas.

Siguiendo con el análisis nos damos cuenta que la columna región parece tener muchos más niveles de los que muestra la función summary() por lo que se decidió utilizar la función levels() para un mejor entendimiento de esta.

levels(avocados$region)

##  [1] "Albany"              "Atlanta"             "BaltimoreWashington"
##  [4] "Boise"               "Boston"              "BuffaloRochester"   
##  [7] "California"          "Charlotte"           "Chicago"            
## [10] "CincinnatiDayton"    "Columbus"            "DallasFtWorth"      
## [13] "Denver"              "Detroit"             "GrandRapids"        
## [16] "GreatLakes"          "HarrisburgScranton"  "HartfordSpringfield"
## [19] "Houston"             "Indianapolis"        "Jacksonville"       
## [22] "LasVegas"            "LosAngeles"          "Louisville"         
## [25] "MiamiFtLauderdale"   "Midsouth"            "Nashville"          
## [28] "NewOrleansMobile"    "NewYork"             "Northeast"          
## [31] "NorthernNewEngland"  "Orlando"             "Philadelphia"       
## [34] "PhoenixTucson"       "Pittsburgh"          "Plains"             
## [37] "Portland"            "RaleighGreensboro"   "RichmondNorfolk"    
## [40] "Roanoke"             "Sacramento"          "SanDiego"           
## [43] "SanFrancisco"        "Seattle"             "SouthCarolina"      
## [46] "SouthCentral"        "Southeast"           "Spokane"            
## [49] "StLouis"             "Syracuse"            "Tampa"              
## [52] "TotalUS"             "West"                "WestTexNewMexico"

Como se puede observar, existen demasiados niveles en la columna como para poder ser manejada de manera efectiva, por lo que consultando el entendimiento sobre la columna de otras personas que también trabajaron la base de datos pero que son de Estados Unidos, nos percatamos de que esta tiene combinadas regiones generales de Estados Unidos, como norte, sur, etc, con ciudades del país, por lo que se dificulta el análisis si estás se mantienen combinadas. Por lo tanto, se decidió separar la base de datos original en varias que tomen en cuenta estas diferencias.

Nombramiento y separación

av1 <- av

av1 <- rename(av1, small_hass = "X4046", large_hass = "X4225", xl_hass = "X4770")

av1_region <- av1 %>%
  filter(region %in% c("California", "West", "SouthCentral", "GreatLakes", "Midsouth", "Southeast", "Northeast", "Plains"))

 # av2_region<- av1 %>%
 #   filter(region %in% c("California", "West", "SouthCentral", "GreatLakes", "Midsouth", "Southeast", "Northeast", "Plains"))

av1_market <- av1 %>%
  filter(!(region %in% c("California", "West", "SouthCentral", "GreatLakes", "Midsouth", "Southeast", "Northeast", "Plains", "TotalUS")))

av1_total <- av1 %>%
  filter(region == "TotalUS")

Análisis de precios según su tipo

Una vez separadas las regiones y ciudades de Estados Unidos, pasamos al análisis del comportamiento de precios. Iniciando por identificar la distribución de los precios por tipo de aguacate, creamos una tabla que grafique la densidad de Kernel dado que un histograma normal sería más complicado de interpretar dada la gran cantidad de precios que se piensan graficar, y al ser un método no paramétrico no supone un riesgo de sesgo al momento de analizar la gráfica.

ggplot(av1_region ,aes(x = AveragePrice, fill = type)) +
  geom_density(alpha = 0.7) +
  geom_vline(xintercept = 1.36, linetype = "dashed") +
scale_x_continuous(breaks = seq(0, 3, 0.3)) +
scale_fill_manual(values = c("#356211","#FFC324")) +
theme_minimal() +
labs(title = "Distribución de los Precios del Aguacate por Tipo",
       x = "Precio Promedio",
       y = "",
       fill = "Tipo de Aguacate") 

Como podemos ver en la gráfica, los aguacates orgánicos tienden a ser más caros que los convencionales como era de esperarse, sin embargo, su comportamiento en la variación de sus precios es diferente, pues mientras que los aguacates orgánicos parecen tener una mayor variación en sus precios, los convencionales se mantienen en su mayoría en el mismo rango. Pero ¿es está una diferencia que afecte a la venta de aguacates y por lo tanto a nuestras predicciones?

Ventas de aguacate según su tipo

Para identificar si las variaciones de precios entre los tipos de aguacates tienen un efecto en sus ventas, graficamos las ventas de aguacate divididas en tipos por los 4 años que hay en la base de datos.

av1_region %>%
ggplot(aes(year, Total.Volume, fill = type)) + 
    geom_bar(width = 0.5, stat = "identity") +
scale_x_continuous(breaks = seq(2015, 2018, 1)) +
scale_y_continuous(labels = label_number(suffix = "B", scale = 1e-8)) +
scale_fill_manual(values = c("#356211", "#ffc324")) +
theme_minimal() +
labs(title = "Ventas de Aguacate",
       x = "Año",
       y = "Volumen de venta en Billones",
       fill = "Tipo")

Como se puede apreciar, la venta de aguacates convencionales supera por mucho las ventas del orgánico, por lo que efectivamente el que el convencional tenga un costo menor supone una ventaja para la venta de aguacates. Esto podría significar que para la predicción efectiva de las ventas de aguacate, lo mejor será dividir o directamente descartar a los aguacates orgánicos, pues podrían suponer outliers innecesarios. Otro hallazgo de la gráfica es que la venta de aguacates no parece tener temporalidad de tipo anual, sin embargo, esto no descarta temporalidad de tipo mensual.

Análisis de temporalidad mensual

Para poder identificar temporalidad, vamos a graficar las ventas totales en todo EU a través de los meses durante los 3 años y 4 meses de los que tenemos datos.

av1_total$Date <- as.Date(av1_total$Date, "%d/%m/%Y")

ggplot(av1_total ,aes(Date, Total.Volume, color = type, group = type)) +
    geom_line(size = 0.8) +
scale_x_date(date_labels = "%b-%y",
               date_breaks =  "3 month",
               limits = c(min(date(av1_total$Date)), max(date(av1_total$Date)))) +
scale_y_continuous(breaks = waiver()) +
scale_color_manual(values = c("#356211", "#ffc324")) +
theme_minimal() +
labs(title = "Volumen de ventas por mes en todo EU",
       x = "Fecha",
       y = "Volumen",
       color = "Tipo") 

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Como claramente se puede apreciar en la gráfica, existe una temporalidad en las ventas de aguacate convencional, el cual se suele vender mucho más en febrero y decae a su punto más bajo en noviembre. Esta temporalidad probablemente se debe al super bowl. Otro punto a rescatar de la gráfica es que los aguacates orgánicos no cuentan con temporalidad, lo cual es otra razón más para separarlos al momento de hacer las predicciones.

Análisis de precios mensuales

Siguiendo con el anterior análisis, ahora graficamos de la misma manera, pero para los precios promedios, buscando identificar correlación con la temporalidad de ventas.

av1_total$Date <- as.Date(av1_total$Date, "%d/%m/%Y")

ggplot(av1_total ,aes(Date, AveragePrice, color = type, group = type)) +
    geom_line(size = 0.8) +
scale_x_date(date_labels = "%b-%y",
               date_breaks =  "3 month",
               limits = c(min(date(av1_total$Date)), max(date(av1_total$Date)))) +
scale_y_continuous(breaks = waiver()) +
scale_color_manual(values = c("#356211", "#ffc324")) +
theme_minimal() +
labs(title = "Precios promedios por mes en todo EU",
       x = "Fecha",
       y = "Precio Promedio",
       color = "Tipo") 

En está gráfica también podemos encontrar temporalidad, sin embargo y a nuestro asombro, es una temporalidad contraría a la esperada, pues en lugar de que los precios sean altos durante el superbowl cuando la demanda es mayor, resulta que son los más bajos del año, lo que si era esperado es que en noviembre, cuando se vende menos, los precios sean los más altos.

Análisis de ventas por región

Siguiendo con el análisis de los datos, ahora pasamos a las ventas de aguacate por región. Además, a este análisis se le va incluir la distribución de estas ventas en los diferentes tipos de “Bags” que vienen en la base de datos, para esto, se transformaron las columnas de “Bags” a filas con sus respectivos números a través de la función pivot_longer()

av1 <- av1 %>%
  pivot_longer(c(Small.Bags, Large.Bags, XLarge.Bags), names_to = "bag_size", values_to = "bag_total")

ggplot(av1_region, aes(reorder(region, Total.Volume), Total.Volume, fill = bag_size)) +
geom_bar(stat = "identity", width = 0.7) +
scale_y_continuous(labels = label_number(suffix = "B", scale = 1e-9)) +
scale_fill_manual(values = c("#356211", "#ffc324", "#e3655b", "#3a3335")) +
coord_flip() +
theme_minimal() +
labs(title = "Ventas de aguacate por región y bag",
       x = "Región",
       y = "Volumen de ventas en Billones",
       fill = "Tipo de Bag")

Como se puede ver, la región con más ventas fue West con más de 3 billones y la que tuvo menos fue la de Plains con poco menos de 1 billón de ventas, más en este aspecto no parece haber un descubrimiento que nos obligue a modificar la base de datos antes de hacer el modelo de predicciones, así como también pudimos descubrir que en realidad el tipo de bolsa o “Bag” realmente parece seguir un patrón normal donde se distribuyen equitativamente entre todas las ventas, por lo que tampoco suponen una anormalidad que se pueda considerar. Sin embargo, para estar seguros ahora graficamos un mapa de volumen de ventas por ciudad/market.

Mapa de volumen de ventas.

library(maps)

## 
## Attaching package: 'maps'

## The following object is masked from 'package:purrr':
## 
##     map

library(mapproj)

states<- map_data("county")
 names(states)[names(states) == "region"] <- "region2"
 names(states)[names(states) == "subregion"] <- "region"
ar <- av1_market



ar$region<- names(ar$region)<-tolower(ar$region)

ar.geo <- merge(states,ar, sort = F, by = "region", all.x = T)
ar.geo<- ar.geo[order(ar.geo$order),]

ggplot(ar.geo, aes(long, lat))+
  geom_polygon(aes(group= group, fill=Total.Volume), color= "white")+
  scale_fill_continuous(
    low = "blue"
      ,high="red",  na.value = "white",
    guide= "colorbar",
    
  )+
  coord_map()

ar2 <- av1_region

ar2$region<- names(ar2$region)<-tolower(ar2$region)

# Get the map data for the US states
ar.geo2 <- map_data("state")

# Merge the map data with the arrests data
ar.geo2 <- merge(ar.geo2, ar2, by.x = "region", all.x = TRUE)

# Create the plot
ggplot(ar.geo2, aes(long, lat)) +
  geom_polygon(aes(group = group, fill = Total.Volume), color = "white") +
  scale_fill_continuous(low = "blue", high = "red", na.value = "blue", guide = "colorbar") +
  coord_map()

Modelos de predicción

Preparación de los datos

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

# modelo de prediccion arima
ts_av <- ts(av$Total.Volume, frequency = 12, start = c(2014, 1), end = c(2100,12))

Modelo de predicción ARIMA

auto.arima(ts_av)

## Series: ts_av 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1      mean
##       0.9778  -0.2784  809099.7
## s.e.  0.0065   0.0300  384307.5
## 
## sigma^2 = 1.404e+11:  log likelihood = -14879.71
## AIC=29767.43   AICc=29767.46   BIC=29787.23

model_av <- arima(ts_av, order = c(1,1,2))

summary(model_av)

## 
## Call:
## arima(x = ts_av, order = c(1, 1, 2))
## 
## Coefficients:
##           ar1      ma1      ma2
##       -0.2218  -0.0609  -0.0845
## s.e.   0.3143   0.3131   0.0926
## 
## sigma^2 estimated as 1.415e+11:  log likelihood = -14869.81,  aic = 29747.61
## 
## Training set error measures:
##                    ME   RMSE      MAE       MPE     MAPE      MASE         ACF1
## Training set 55.03804 375985 119058.9 -12.62595 25.31166 0.9633199 -0.004317561

checkresiduals(model_av)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,2)
## Q* = 130.92, df = 21, p-value < 2.2e-16
## 
## Model df: 3.   Total lags used: 24

forecast_av <- forecast(model_av, h = 12)

plot(forecast_av, main = "Predicción de volumen total hasta 2100", xlab = "Año", ylab = "Total.Volume", xlim = c(2020, 2100), ylim = c(0, 200000), type = "l", col = "blue")

#Predicción del modelo ETS

#prediccion del modelo ets 

ts_av2 <- ts(av$Total.Volume, frequency = 12, start = c(2020, 1), end = c(2030, 12))
 
ets_model <- ets(ts_av2)

forecast_av <- forecast(ets_model, h = 120)

plot(forecast_av, main = "Prediccion del volumen total hasta 2030", xlab = "Año", ylab = "Total.Volume")

Se utlizaron dos modelos de prediccion uno con arima y otro con la tecnica de suavizado exponencial etc ya que creemos que estos dos modelos son los mas optimos para predecir una serie temporal con tendencias para la variable del volumen total y nos ayudara a poder planear a futuro