RETO:
El reto se plantea bajo la siguiente pregunta: A partir de los datos abiertos del portal MEDATA, ¿pueden las aerolíneas ayudarles a los hoteles a predecir su demanda?
Tras la búsqueda de información en el portal MEDATA, se encontraron las siguientes bases de datos que pueden ayudar a resolver el reto planteado:
- BD con datos de llegada mensual de pasajeros al aeropuerto con origen internacional.
- BD con datos de llegada mensual de pasajeros al aeropuerto con origen nacional.
- BD con datos porcentuales de ocupación hotelera mensual de Medellín.
Con la anterior información, se sugiere transformar la pregunta del reto, sin alterar su escencia: con base en los datos de llegada de personas vía aerea a la ciudad, ¿es posible predecir el porcentaje de ocupación hotelera de la ciudad?
Cargamos las librerías necesarias para el reto:
Cargamos las bases de datos identificadas en el portal MEDATA:
# - Cargamos archivo con información de vuelos internacionales
vuelos_i <- read.csv("~/Dropbox/base_vuelos_hoteles/llegada_mensual_pasajeros_aeropuerto_de_origen_internacional.csv", header = TRUE, sep = ';')
# - Cargamos archivo con información de vuelos nacionales
vuelos_n <- read.csv("~/Dropbox/base_vuelos_hoteles/llegada_pasajeros_mensual_por_aeropuerto_de_origen_nacional.csv", header = TRUE, sep = ';')
# - Cargamos archivo con informaciòn de ocupación hotelera de la ciudad
ocupacion_h <- read.csv("~/Dropbox/base_vuelos_hoteles/porcentaje_ocupacion_hotelera_mensual_de_medellin.csv", header = TRUE, sep = ';')- La base de datos de vuelos internacionales tiene 4861 filas y 5 columnas, reporta datos desde el 200701 hasta el 201906 y contiene en total las siguientes variables:
## [1] "periodo" "nombre" "codigo" "indicador" "valor"- La base de datos de vuelos nacionales tiene 10549 filas y 4 columnas, reporta datos desde el 200701 hasta el 201906 y contiene en total las siguientes variables:
## [1] "periodo" "codigo" "indice" "valor"- La base de datos de ocupación hotelera tiene 153 filas y 3 columnas, reporta datos desde el 200701 hasta el 201909 y contiene en total las siguientes variables:
## [1] "periodo" "indice" "valor"Ahora, debemos de crear un nuevo dataframe que nos agrupe toda la información de los vuelos nacionales e internacionales con la ocupación hotelera teniendo en centa solo las variables relevantes para el ejercicio:
# - Realizamos un subset para seleccionar solo las variables que nos interesan para el reto:
vuelos_i <- subset(vuelos_i, select = c("periodo", "valor"))
vuelos_n <- subset(vuelos_n, select = c("periodo", "valor"))
ocupacion_h <- subset(ocupacion_h, select = c("periodo", "valor"))
# - Realizamos una agrupación para consolidar los vuelos nacionales e internacionales
consolidado_vi <- vuelos_i %>% group_by(periodo) %>% summarise(total=sum(valor))## `summarise()` ungrouping output (override with `.groups` argument)## `summarise()` ungrouping output (override with `.groups` argument)# - Renombramos variables
names(consolidado_vi) <- c("periodo", "total_vi")
names(consolidado_vn) <- c("periodo", "total_vn")
names(ocupacion_h) <- c("periodo", "ocupacion")
# - Filtramos en todas las BD, el máximo periodo con información completa en las tres BD.
consolidado_vi <- subset(consolidado_vi, subset = (periodo < 201904))
consolidado_vn <- subset(consolidado_vn, subset = (periodo < 201904))
ocupacion_h<-subset(ocupacion_h, subset = (periodo < 201904))
# utilizamos la función Merge para unir las bases de datos
ocupacion_vuelos <- merge(x = consolidado_vi, y = consolidado_vn, by = "periodo")
datos_an <- merge(x = ocupacion_vuelos, y = ocupacion_h, by = "periodo")Realizamos una visualización para identificar las correlaciones
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(datos_an,lower.panel = panel.smooth,upper.panel = panel.cor)
Teniendo en cuenta que las correlaciones entre las variables periodo, total_vi, y total_vn tienen una alta correlación, seleccionaremos solo la variable total_vn para correr un modelo de regresión lineal teniendo como variable dependiente la ocupación:
Adicionalmente, es importante indicar que se verificó una tabla de correlacionescon transformación de variables a escala logarítmica, pero lamentablemente no se evidenció una mejor dispersión o correlación de los datos.
Por lo anterior, se decide trabajar con las variables sin realizar las transformaciones logaritmicas y por tanto tampoco realizamos escalamiento de datos.
Transformando la variable Periodo
Tras analizar la variable periodo, se considera pertinente transformarla para tener el año y el mes en variables separadas y tratarlas como variables categóricas.
# - Separamos la variable periodo en Año y Mes
datos_an <- separate(datos_an, periodo, c("anio", "mes"),
sep = 4, convert = T)
# - Verificamos con la función names
names(datos_an)## [1] "anio" "mes" "total_vi" "total_vn" "ocupacion"Para poder incluir estas variables categóricas dentro del modelo, debemos transformarlas en variables de clase factor. Procedemos:
Calibración de un modelo lineal
set.seed(0574)
p_vl <- 0.3
n <- dim(datos_an)[1]
ix_vl<-sample(n,round(0.3*n),replace = FALSE)
datos_tr<-datos_an[-ix_vl,]
datos_vl<-datos_an[ix_vl,]Ajuste del modelo Lineal
##
## Call:
## lm(formula = datos_tr$ocupacion ~ datos_tr$total_vn + datos_tr$anio +
## datos_tr$mes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9069 -1.7500 -0.1985 1.4211 9.9569
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.685e+01 2.899e+00 16.161 < 2e-16 ***
## datos_tr$total_vn 9.405e-05 2.430e-05 3.870 0.000225 ***
## datos_tr$anio2008 -5.425e+00 1.692e+00 -3.207 0.001946 **
## datos_tr$anio2009 -1.217e+01 1.705e+00 -7.134 4.37e-10 ***
## datos_tr$anio2010 -1.307e+01 2.028e+00 -6.446 8.74e-09 ***
## datos_tr$anio2011 -1.303e+01 2.077e+00 -6.276 1.81e-08 ***
## datos_tr$anio2012 -1.402e+01 2.833e+00 -4.946 4.25e-06 ***
## datos_tr$anio2013 -1.597e+01 4.260e+00 -3.748 0.000340 ***
## datos_tr$anio2014 -1.245e+01 3.988e+00 -3.123 0.002515 **
## datos_tr$anio2015 -1.204e+01 4.383e+00 -2.747 0.007459 **
## datos_tr$anio2016 -1.026e+01 4.778e+00 -2.147 0.034939 *
## datos_tr$anio2017 -1.401e+01 4.466e+00 -3.136 0.002412 **
## datos_tr$anio2018 -1.761e+01 5.182e+00 -3.397 0.001075 **
## datos_tr$anio2019 -1.201e+01 5.766e+00 -2.083 0.040544 *
## datos_tr$mes2 7.076e-01 1.756e+00 0.403 0.688147
## datos_tr$mes3 -5.345e-01 1.707e+00 -0.313 0.755049
## datos_tr$mes4 1.303e+00 1.881e+00 0.692 0.490695
## datos_tr$mes5 6.887e-01 1.722e+00 0.400 0.690243
## datos_tr$mes6 -4.441e-01 1.641e+00 -0.271 0.787430
## datos_tr$mes7 2.861e+00 1.703e+00 1.680 0.096985 .
## datos_tr$mes8 8.683e+00 1.794e+00 4.839 6.45e-06 ***
## datos_tr$mes9 4.662e+00 1.829e+00 2.549 0.012754 *
## datos_tr$mes10 3.684e+00 1.727e+00 2.134 0.036015 *
## datos_tr$mes11 6.560e+00 1.772e+00 3.702 0.000397 ***
## datos_tr$mes12 1.341e+00 1.865e+00 0.719 0.474316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.105 on 78 degrees of freedom
## Multiple R-squared: 0.8609, Adjusted R-squared: 0.8181
## F-statistic: 20.11 on 24 and 78 DF, p-value: < 2.2e-16
Los resultados de este ejercicio nos muestra un mejor resultado ya que el R cuadrado alcanza un 0.86. por tanto, se considera que con los datos disponibles es el porcentaje más alto que se puede alcanzar con este modelo.
Gráficos diagnósticos:
Residuales versus valores predichos:
En el gráfico se pude ser que no hay patrones exagerados aunque hay observaciones sospechosas como la 32, la 6 y el 85.
Gráfico cuantil-cuantil de los residuales estandarizados vs los de una distribución normal:
Vemos que al final en la cola hay unos datos que se despegan de la recta identidad por tanto puede que ahí no se esten cumpliendo los supuestos de normalidad
Residuales estandarizados versus valores predichos:
Este gráfico confirma los residuales sospechosos que se identificaron anteriormente y vemos que los datos están por debajo de 2 lo cual es bueno.
Distancia de Cook de cada observación:
Aquí podemos ver las observaciones con distancia de Cook altas.
Residuales estandarizados versus leverage:
Distancia de Cook vs Leverage:
Aquí vemos las observaciones que están influyendo mucho en el modelo y aumentan los errores de predicción.