library(dplyr)
library(ggplot2)
library(prophet)
library(lubridate)
El dataset consiste en información real de 8 meses de visitantes de un local de ropa localizado en un shopping de gran tamaño
shopping <- read.csv('local-shopping_prophet.csv') %>%
rename(.,visitantes=Total.de.Visitantes) %>%
mutate(Dia = ymd(Dia))
glimpse(shopping)
Observations: 228
Variables: 3
$ Dia [3m[38;5;246m<date>[39m[23m 2017-09-14, 2017-09-15, 2017-09-16, 2017-09-17, 2017-09-18, 2017-09-19, 2017-09-20, 2017-09-21, 2017-09-22…
$ visitantes [3m[38;5;246m<int>[39m[23m 589, 696, 1034, 940, 540, 526, 312, 553, 639, 793, 745, 8, 333, 568, 352, 393, 497, 578, 563, 520, 413, 8, …
$ Clima [3m[38;5;246m<fct>[39m[23m Parcialmente Nublado, Soleado, Soleado, Parcialmente Nublado, Soleado, Soleado, Soleado, Nublado, Soleado, …
Se observa que la serie de tiempo presenta un comportamiento cíclico con algunos eventos atípicos:
ggplot(shopping, aes(Dia, visitantes)) + geom_line() + theme_bw() + labs(title='Visitantes por dia')
ggplot(shopping, aes(Dia, visitantes)) + geom_line() + geom_point(color='forestgreen') + theme_bw() + labs(title='Visitantes por dia')
Una primera aproximación puede ser tratar de utilizar técnicas de suavizado.
En este caso estamos utilizando loess como una herramienta gráfica con el comando geom_smooth()
de ggplot.
ggplot(shopping, aes(Dia, visitantes)) + geom_point(color='forestgreen') + geom_smooth() + theme_bw() + labs(title='Visitantes por dia: Suavizado')
Si bien está suavizando, el resultado es bastante pobre ya que el modelo no está logrando representar bien la variabilidad de los datos.
Ahora realizamos 4 modelos de LOESS modificando la ventana (span) para cada uno de ellos.
loess=stats::loess(visitantes~as.numeric(Dia), data = shopping, na.action = 'na.exclude', model = T, span=0.01)
shopping['loess']=predict(loess,shopping)
loess_1 = ggplot(shopping, aes(Dia,visitantes)) + geom_point() + geom_line(aes(y=loess), color='firebrick', size=1) +
labs(title= "LOESS span:0.01") + theme_bw()
loess=stats::loess(visitantes~as.numeric(Dia), data = shopping, na.action = 'na.exclude', model = T, span=0.25)
shopping['loess']=predict(loess,shopping)
loess_2 = ggplot(shopping, aes(Dia,visitantes)) + geom_point() + geom_line(aes(y=loess), color='forestgreen', size=1) +
labs(title= "LOESS span:0.25") + theme_bw()
loess=stats::loess(visitantes~as.numeric(Dia), data = shopping, na.action = 'na.exclude', model = T, span=0.5)
shopping['loess']=predict(loess,shopping)
loess_3 = ggplot(shopping, aes(Dia,visitantes)) + geom_point() + geom_line(aes(y=loess), color='steelblue', size=1) +
labs(title= "LOESS span:0.50") + theme_bw()
loess=stats::loess(visitantes~as.numeric(Dia), data = shopping, na.action = 'na.exclude', model = T, span=0.75)
shopping['loess']=predict(loess,shopping)
loess_4 = ggplot(shopping, aes(Dia,visitantes)) + geom_point() + geom_line(aes(y=loess), color='purple', size=1) +
labs(title= "LOESS span:0.75") + theme_bw()
cowplot::plot_grid(loess_1, loess_2, loess_3,loess_4)
Se nota que para un span muy pequeño nuestro modelo va a tener un problema de OVERFITTING: no va a poder generalizar. Mientras que para spans grandes tenemos un problema de UNDERFITTING: nuestro modelo no logra captar bien la variabilidad del fenómeno.
Vamos a trabajar sobre la implementación de Prophet en R.
Prophet requiere que le pasemos el dataset con:
ds: la variable temporal
y: la variable a predecir
# Eliminamos observaciones con menos de 250 visitantes (cuestion de negocio)
shopping[shopping['visitantes']<250,'visitantes'] = NA
# Creamos el dataset
prophet_df = shopping %>% select(Dia, visitantes) %>% rename(., ds=Dia, y=visitantes)
Recordemos que el modelo subyacente es el propuesto por Harvey & Peters (1990):
\(y(t)=g(t)+s(t)+h(t)+\epsilon_t\)
Hay muchos parametros para tener en cuenta. Veremos algunos:
df: dataframe
growth: tipo de tendencia: lineal o logistica
yearly.seasonality: hay estacionalidad anual?
yearly.seasonality: hay estacionalidad diaria?
holidays: dataframe con fechas de vacaciones/eventos especiales
La funcion prophet
crea el modelo, podemos pasarle o no el dataframe.
La funcion prophet.fit
aplica un modelo creado a un dataframe
# Llamamos solo al modelo
prophet_base=prophet()
# Le pasamos el dataset
prophet_base = fit.prophet(m = prophet_base, prophet_df)
Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Initial log joint probability = -4.68957
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Notemos que el modelo automaticamente deshabilita la estacionalidad anual y diaria.
Llamando a plot
obtenemos el valor predicho del modelo y el valor original.
Es dificil de notar pero el modelo realiza predicciones aun para los dias en los cuales no hay datos.
plot(prophet_base,fcst=predict(prophet_base, prophet_df)) +theme_bw()
La funcion prophet_plot_components
nos devuelve los efectos de los componentes en nuestra variable a predecir
prophet_plot_components(prophet_base, fcst=predict(prophet_base, prophet_df)) +theme_bw()
NULL
En este caso tenemos la tendencia de los visitantes y la estacionalidad semanal. Esta última parece tener bastante sentido ya que nos muestra que hay mayor cantidad de visitantes durante el fin de semana.
La funcion cross_validation
permite realizar pronosticos realizando un esquema de cross-validation temporal y, a partir de ellos, obtener ciertas metricas de performance.
horizon: horizonte del pronostico. Cuánto tiempo deseo predecir
period: periodo entre fechas de analisis. Cuánto voy corriendo mi ventana
initial: periodo inicial de entrenamiento
cv_base = cross_validation(prophet_base, initial = 45, period = 7, horizon = 15, units = 'days')
Making 24 forecasts with cutoffs between 2017-11-04 and 2018-04-14
Initial log joint probability = -3.06265
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -2.89179
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.42994
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -2.99569
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.58187
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.40797
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -6.14203
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -6.67306
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.38137
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.21596
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.43862
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.36873
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.43537
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.5117
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.60968
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.69348
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.77909
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.8066
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -4.02179
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -4.02246
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -3.99051
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -4.04036
Optimization terminated normally:
Convergence detected: relative gradient magnitude is below tolerance
Initial log joint probability = -4.27091
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Initial log joint probability = -4.33551
Optimization terminated normally:
Convergence detected: relative gradient magnitude is below tolerance
cv_base
La funcion performance_metrics
computa varias metricas de performance a partir de un dataframe de cross validation de prophet
performance_metrics(cv_base, rolling_window = 0.5)
La funcion add_seasonality
nos permite agregar nuevas estacionalidades. Las estacionalidades se modelan utilizando series de Fourier que nosotros debemos definir.
Definimos:
m: modelo
name: nombre de la estacionalidad
period: cantidad de dias del periodo
fourier.order: orden de la serie de fourier para modelar la estacionalidad
# Llamamos solo al modelo
prophet_mensual=prophet()
# Agregamos la estacionalidad mensual
prophet_mensual=add_seasonality(prophet_mensual, name='monthly', period=365/12, fourier.order = 4)
# Le pasamos el dataset
prophet_mensual = fit.prophet(m = prophet_mensual, prophet_df)
Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Initial log joint probability = -4.68957
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Notemos que sigue deshabilitando la estacionalidad anual y diaria.
plot(prophet_mensual,fcst=predict(prophet_mensual, prophet_df)) +theme_bw()
La funcion prophet_plot_components
nos devuelve los efectos de los componentes en nuestra variable a predecir
prophet_plot_components(prophet_mensual, fcst=predict(prophet_mensual, prophet_df)) +theme_bw()
NULL
Notemos que la tendencia y estacionalidad semanal se mantiene y ahora se agrega la tendencia mensual que habíamos incorporado. Nos indica que el ciclo mensual se caracteriza por un crecimiento en la mitad de mes.
La funcion performance_metrics
computa varias metricas de performance a partir de un dataframe de cross validation de prophet
performance_metrics(cv_mensual, rolling_window = 0.5)
Como último paso vamos a agregar las ventas de navidad y ciertos dias de promociones como eventos especiales
Creamos el dataframe de eventos con: nombre del evento, fechas y una “ventana” para definir si el evento se estira a ciertos dias.
# Navidad
christmas = data.frame(holiday= 'christmas',
ds=ymd(c('2017-12-16','2017-12-17','2017-12-18',
'2017-12-19','2017-12-20','2017-12-21',
'2017-12-22','2017-12-23')),
lower_window= 0,
upper_window= 0)
# Promociones
big_sales = data.frame(
holiday= 'big_sales',
ds= ymd(c('2017-09-16','2017-10-08','2017-10-14',
'2017-11-20','2017-12-03','2017-12-30')),
lower_window= 0,
upper_window= 0)
holidays= bind_rows(christmas, big_sales)
glimpse(holidays)
Observations: 14
Variables: 4
$ holiday [3m[38;5;246m<chr>[39m[23m "christmas", "christmas", "christmas", "christmas", "christmas", "christmas", "christmas", "christmas", "…
$ ds [3m[38;5;246m<date>[39m[23m 2017-12-16, 2017-12-17, 2017-12-18, 2017-12-19, 2017-12-20, 2017-12-21, 2017-12-22, 2017-12-23, 2017-09-…
$ lower_window [3m[38;5;246m<dbl>[39m[23m 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
$ upper_window [3m[38;5;246m<dbl>[39m[23m 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
# Llamamos al modelo con el dataset de eventos
prophet_full=prophet(holidays = holidays)
# Agregamos la estacionalidad mensual
prophet_full=add_seasonality(prophet_full, name='monthly', period=30.5, fourier.order = 4)
# Le pasamos el dataset
prophet_full = fit.prophet(m = prophet_full, prophet_df)
Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Initial log joint probability = -4.68957
Optimization terminated normally:
Convergence detected: absolute parameter change was below tolerance
Fijense que el modelo automaticamente deshabilita la estacionalidad anual y diaria.
plot(prophet_full,fcst=predict(prophet_full, prophet_df)) +theme_bw()
La funcion prophet_plot_components
nos devuelve los efectos de los componentes en nuestra variable a predecir
prophet_plot_components(prophet_full, fcst=predict(prophet_full, prophet_df)) +theme_bw()
NULL
La tendencia y estacionalidad semanal se mantienen aproximadamente igual. La estacionalidad mensual cambia bastante y se observa que existen tres picos en el ciclo mensual. Por su parte, los eventos se modelan como pequeños saltos o picos.
performance_metrics(cv_full, rolling_window = 0.5)
Prophet también nos permite realizar gráficos interactivos que suelen ser muy útiles para presentar los resultados.
dyplot.prophet(prophet_full, fcst=predict(prophet_full, prophet_df))