Este trabajo tiene como motivo principal el mostrar ciertos recursos visuales y de análisis de datos para validar capacidades en el lenguaje de programación R. Para esta secuencia de trabajos, utilize los datos de las exportaciones realizadas por Chile entre 2010 y 2020, datos entregados por la subsecretaría de transporte como parte de un proceso de selección laboral del que forme parte entre mayo y junio del presente año.
Los datos originales fueron facilitados mediante claves de acceso a instancia de base de datos Postgres. Mediante conexión a esta información via API y luego limpieza y consolidación de los datos, procedí a enriquecer los datos mediante conexiones también via API al banco mundial para los datos de PIB y climáticos (mm de lluvia para tener algún factor demográfico para explorar), para lo cual fue necesario, realizar labores de traducción apoyándome de data externa. Posteriormente, alimente los datos con información de geolocalización para realizar proyecciones visuales con los datos de latitud y longitud adquiridos, en primaria instancia, con los vectores regionales de la biblioteca nacional, pero al confirmar cierta imprecisión en los datos, procedí a utilizar los geo-datos del INE.
El trabajo consistía en realizar proyecciones, algoritmos de predicción y de clustering, además de trabajo de exploración y presentación.
Posterior al período de entrega de ciertos requerimientos asociados a estos datos para fines de la selección laboral (4 días en primera instancia y luego 3 días más), continué explorando los datos por iniciativa propia, desarrollando variadas visualizaciones y algoritmos como parte de mi aprendizaje y luego para enseñar para futuras postulaciones a otros cargos, y en eso consiste este documento.
library(data.table)
library(ggthemes)
library(tidyverse)
library(ggthemes)
library(ggrepel)
library(gridExtra)
library(skimr)
library(zoo)
library(ggpmisc)
library(kableExtra)
library(ggpubr)
library(stringr)
library(forcats)
library(treemapify)
library(zoo)
library(univariateML)
library(wesanderson)#Lectura inicial de los datos. Todas las columnas del tipo “character” son convertidas al tipo factor.
options(scipen= 999)
datos <- fread("archivo1_.csv")
if ("V1" %in% names(datos)) {datos[, c('V1')] <- NULL}
setDT(datos)
cambioCol <- names(Filter(is.character, datos))
datos[, (cambioCol):= lapply(.SD, as.factor), .SDcols = cambioCol]
glimpse(datos)## Rows: 2,668,403
## Columns: 27
## $ año <int> 2011, 2011, 2011, 2011, 2011, 20...
## $ mes <int> 10, 10, 10, 10, 10, 10, 10, 10, ...
## $ cod_tipo_carga_operacionexpo <fct> F, R, R, R, R, F, R, R, R, F, F,...
## $ item_sa_operacionexpo <int> 3041942, 3021221, 3041942, 60319...
## $ valor <dbl> 18311.25, 15667.26, 6023.99, 215...
## $ peso <dbl> 2695, 3272, 1029, 1403, 11100, 6...
## $ glosa_regionorigen <fct> REGIÓN DE LOS LAGOS, REGIÓN DE L...
## $ nombre_tipo_operacion <fct> EXPORTACIÓN NORMAL, EXPORTACIÓN ...
## $ nombre_aduana <fct> Metropolitana, Metropolitana, Me...
## $ glosa_viatransporte <fct> AÉREO, AÉREO, AÉREO, AÉREO, AÉRE...
## $ nombre_puerto_embarque <fct> AEROP. A.M. BENITEZ, AEROP. A.M....
## $ tipo_puerto_embarque <fct> Aeropuerto, Aeropuerto, Aeropuer...
## $ nombre_puerto_desembarque <fct> OTROS PUERTOS DE PERÚ NO ESPECIF...
## $ pais_puerto_desembarque <fct> Perú, Perú, Perú, Perú, Argentin...
## $ zona_geografica_puerto_desembarque <fct> América del Sur, América del Sur...
## $ nombre_pais <fct> Perú, Perú, Perú, Perú, Argentin...
## $ continente_pais <fct> América, América, América, Améri...
## $ ingles <fct> Peru, Peru, Peru, Peru, Argentin...
## $ long <dbl> 112.74721, 112.74721, 112.74721,...
## $ lat <dbl> 16.65361, 16.65361, 16.65361, 16...
## $ pib <dbl> 6453.561, 6453.561, 6453.561, 64...
## $ lluvia <int> 1738, 1738, 1738, 1738, 591, 591...
## $ lval <dbl> 9, 9, 8, 9, 10, 8, 8, 5, 9, 8, 1...
## $ lpes <dbl> 7, 8, 6, 7, 9, 6, 7, 4, 6, 7, 9,...
## $ fecha <fct> oct. 2011, oct. 2011, oct. 2011,...
## $ estac <fct> pri, pri, pri, pri, pri, pri, pr...
## $ columna <fct> sur, sur, sur, sur, sur, sur, su...
Las variables cualitativas relevantes para este estudio son: - cod_tipo_carga_operacionexpo: el tipo de la carga (exportación) - glosa_regionorigen: la región de origen de la exportación - nombre_puerto_embarque: El puerto donde embarca la exportación - nombre_puerto_desembarque - pais_puerto_desembarque - zona_geografica_puerto_desembarque - nombre_pais: del país receptor de la mercancía - continente_pais: el continente de dicho país
Se explorará en otro apartado la variable estación que cree para poder estudiar el efecto de las estaciones de los países en el número de exportaciones (no estaba incluida en la data original). Ingles es el nombre en inglés del país de destino. La cree para poder hacer match con geo-datos para poder realizar trabajos de geo-exploración. En general la data externa viene en idioma inglés por tanto es bastante conveniente. Pero del total de 218 solo se lograron traducir con cruce de bases de datos gratuitas 159. Hay nombres que no tienen equivalente a nombre de país oficial. Con esta data se pueden luego mapear -que es un trabaje que agregare mas adelante- las exportaciones y trabajar con distancias entre otros. De esta traducción también luego generé las columnas PIB y mm lluvia con conexión a API del bco mundial. Posteriormente, esta data será refinada utilizando el paquete de R rnaturalearth, que contiene datos relacionados a población por censo reciente, categorías de desarrollo económico y de poder adquisitivo.
Uno podría pensar que esta información sería mejor conseguirla directamente del bco mundial, pero tuve malas experiencias por inconsistencia en la disponibilidad de la data via api(a veces se puede descargar otras veces no hay retorno) entre otros detalles.
Otras formas de hacer una data mas robusta para generar diversos análisis incluye la data de exportación mundial que esta disponible creo mediante la ONU (se que está disponible, en el momento de escribir esto no recuerdo exactamente quien la tiene disponible). Esta data incluye todas las exportaciones para un importante número de años, con detalle del producto entre otros.
Otra data que exploré es la de economía regional, que publica el INE, pero tiene bastante data censurada que dificultan el análisis. De todas formas se pueden hacer proxies con los campos disponibles.
Hay muchísimas otras formas de enriquecer los datos, como por ejemplo, crecimiento económico de cada nación en los últimos x meses, proyecciones económicas para estos. Obviamente conociendo mayor detalle del tipo de mercancía exacta uno podría generar estudios bastante sofisticados.
Otra alternativa que exploraré al respecto es la reconversión de la data en datos de panel, donde se agrupan los registros mensuales de cada país y se construye indicadores como valor medio, desviacion, coeficiente de variación, cuartiles de valor, de peso, e indicadores para cada una de las otras variables (valor promedio de las exportaciones de ese país a cada una de las regiones, de las aduanas, etc.) Es decir, pasar de una data de 2.6mm de filas por 28 columnas a una de poco mas de 4000 filas por miles de columnas, donde cada fila sea lo que cierto país exportó en un mes determinado, y cada columnas toda la serie de indicadores señalados. Esto permitiría trabajar la data con toda la gama de características de los datos de panel y también como series de tiempo.
Se re-nivelan las variables nombre_aduana y glosa_regionorigen para que este ordenadas de norte a sur, de modo que cuando se desplieguen en los gráficos mantengan un orden lógico-reconocible.
unicos <- unique(datos[, .(glosa_regionorigen)])[[1]]
regiones <- ifelse(is.na(word(unicos, 4)),
ifelse(is.na(word(unicos, 3)),
word(unicos, 2),
word(unicos, 3)),
word(unicos, 4))
regiones[3] <- "MAGALLANES"
a <- data.frame(unicos, regiones)
datos <- left_join(datos, a, by = c("glosa_regionorigen" = "unicos"))
datos$regiones <- str_to_title(datos$regiones)
datos$regiones <- factor(datos$regiones, levels = rev(str_to_title(c("NACIONALIZADAS", "PARINACOTA", "TARAPACÁ" ,"ATACAMA", "ANTOFAGASTA","COQUIMBO", "VALPARAÍSO", "METROPOLITANA", "O'HIGGINS", "ÑUBLE", "MAULE", "RÍOS", "BIOBIO","ARAUCANÍA", "LAGOS", "AYSÉN", "MAGALLANES"))))Se trabajara sobre una muestrs de 100.000 registros para reducir los costos computacionales.
Los valores lat y long corresponden al resultado del cruce de los datos originales con los datos geo-especiales de fuentes de arcgis. Cada país tiene una latitud y una longitud, por tanto, la razón que explica que sea mas frecuente que la variable nombre_país es que hay diferencias en los decimales no relevantes para efectos de visualización de los mapas. De todas formas, mas adelante se descarta esta variable y se trabaja con los datos del package de R rnaturalearth, que facilita el cálculo de centroides (para mapear con mayor precisión), esta bastante completa y es mas sencilla de invocar, además de que tiene varios indicadores de países como población censal, indicadores de desarrollo económico y otros, que son de utilidad.
Ingles es menor a nombre_país dado que las bases de dato adquiridas para traducción no tenían el término equivalente para la totalidad de los países de la data de la subsecretaría, principalmente por términos como “Otros”, “Territorio Frances en América” etc, que no son nombres totalmente formales.
Lluvia proviene de los datos del banco mundial, y estaban algo incompletos.
lval y lpes son los valores logarítmicos de valor y peso que recalcule sobre los datos originales respectivos, truncando los decimales para que hayan unos pocos grupos (20), que permitirá luego probar algoritmos de clasificación. Este trabajo esta armado y lo subiré con el tiempo.
Aduana y región de origen tienen poco menos de 20 regisotros cada una, y el resto de variables tienen baja cardinalidad.
Las mercaderías son principalmente de la categoría F y R, domina el tipo de puerto de embarque marítimo. Notar que el alto de la barra corresponde al número de veces que aparece cada categoría/tipo de puerto de embarque en una muestra de 100.000 registros o eventos (en particular, exportaciones). Esto representa un 4% de la totalidad de los registros (+2.5mm)
Respecto al continente del país de destino, relacionado a la zona geográfica de desembarcación, queda claro que es América el principal destino de las exportaciones, seguido por Europa y Asia, en términos del número de exportaciones realizadas a dichos destinos.
Hay una gran fortaleza en la zona central respecto al movimiento de mercancias en sus aduanas. Lógicamente, los principales puertos marítimos dominan en este ámbito, en particular los puertos de Valparaíso, San Antonio y Talcahuano. La región metropolitana, probablemente por vía aérea (esto será explorado mas adelante) y Los Andes por vía terrestre. También hay una clara correlación entre el número de habitantes de la región y el número de habitantes, relación que también sera explorada posteriormente, dado que esta es solo un estudio de las variables individuales (con la particularidad, de que pareciera que la región de Ohiggins genera más por habitante que el resto de las comunas, este ratio también sera explorado) Puertos de Embarque
Claramente los principales puertos de embarque se corresponden con el dinamismo de las aduanas, visto en el gráfico anterior. Queda patente la relación 1:1 entre aduana y puerto de embarque (es decir, cada aduana representa principalmente 1 puerto de embarque; San Antonio-San Antonio, Aereop. Benitez - Metropolitana, etc.)
Puertos de desembarque
En la data aparecen muchos registros de nombre de puertos de desembarque del tipo “Otros puertos no Especificados”. ¿Que se puede saber de estas exportaciones?
muestra[str_detect(muestra$nombre_puerto_desembarque, "NO ESPECIFICADOS", negate= FALSE),][, .(.N), .(nombre_puerto_desembarque)] %>%
ggplot() + geom_histogram(aes(x=fct_reorder(tolower(nombre_puerto_desembarque), N), y = N,
fill = N), stat = 'identity') +
theme_minimal() +
scale_fill_gradient2_tableau() +
guides(fill = guide_legend(title = "", title.position = "left")) +
labs(title = "Puertos de desembarque que no están especificados",
x = "", y = "Frecuencia") +
coord_flip() +
scale_x_discrete(label = function(x) stringr::str_sub(x, -27)) +
theme(plot.title = element_text(face = "bold"),
legend.position = c(.6, .55),
legend.background = element_rect(fill = "white", size = .5, linetype = "solid", inherit.blank = TRUE))Existen tanto continentes como países relacionados a los datos no especificados. Y de las que si estan especificadas?
muestra[str_detect(muestra$nombre_puerto_desembarque, "NO ESPECIFICADOS", negate= TRUE),][, .(.N), .(nombre_puerto_desembarque)] %>%
ggplot() + geom_histogram(aes(x=fct_reorder(str_to_title(nombre_puerto_desembarque), N),
y = N, fill = N),
stat = 'identity') +
theme_minimal() +
scale_fill_gradient2_tableau() +
labs(title = "Puertos de desembarque que si están especificados",
x = "", y = "Frecuencia") +
coord_flip() +
scale_x_discrete(label = function(x) stringr::str_sub(x, -16)) +
theme(plot.title = element_text(face = "bold"),
legend.position = c(.5, .8),
legend.background = element_rect(fill = "white", size = .5, linetype = "solid", inherit.blank = TRUE)) +
guides(fill = guide_legend(title = "", title.position = "left"))## [1] 0.44389
Sería una variable interesante de explorar debido a su mayor grado de granularidad, pero es necesario limpiar un poco esta variable para poder explorarla con mayor claridad:
muestra$puerto_desembarque <- ifelse(grepl("NO ESPECIFICADOS", as.character(muestra$nombre_puerto_desembarque)),
str_remove(muestra$nombre_puerto_desembarque, "NO ESPECIFICADOS") %>%
str_remove(., "OTROS PUERTOS") %>%
str_remove(., " DE") %>%
trimws() %>%
paste0(.,"_"), as.character(muestra$nombre_puerto_desembarque))Una vez limpiada la variable, se puede estudiar con mayor claridad. Se procede a graficar la variable para los datos no y si especificados, pero dado que son más de 210, y de que existe un grupo de puertos importantes con apenas registros, se procede luego a filtrar solo aquellos que hayan realizado más de 200 exportaciones (recordar que se esta trabajando con una muestra de 100.000 datos para reducir el costo computacional).
Notar que en el próximo gráfico hay ciertas variables con terminación "_". Esta terminación la cree para que se pueda distinguir aquellos puertos de desembarque que fueron modificados. La idea es que se pueda indentificar aquellos puertos de desembarque que no estaban especificados en la data original, de los que sí.
El principal puerto es de América Central, y le siguen varios de América del Sur y de Estados Unidos, además detacando el puerto de Roterdam, uno asiático no especificado y China en menor medida. Subconjunto de solo los puertos de desembarque con más de 200 exportaciones.
El puerto no especificado de América se refiere concretamente a:
unique(muestra[puerto_desembarque == "AMÉRICA_",
.(zona_desembarque = zona_geografica_puerto_desembarque, pais_puerto_desembarque, .N ),
.(puerto_desembarque)])## puerto_desembarque zona_desembarque pais_puerto_desembarque N
## 1: AMÉRICA_ América Central Antillas Holandesas 13320
Y concretamente, los puertos no especificados de Asia se refieren a:
unique(muestra[puerto_desembarque == "ASIÁTICOS_",
.(zona_desembarque = zona_geografica_puerto_desembarque, pais_puerto_desembarque, .N ),
.(puerto_desembarque)])## puerto_desembarque zona_desembarque pais_puerto_desembarque N
## 1: ASIÁTICOS_ Asia Singapur 3280
Lo siguiente a explorar son los paises de destino de las exportaciones agrupados según la zona geográfica del desembarque. (La variable zona geográfica del desembarque será la usada porque tiene mayor desgloce que la variable continente país (y por tanto más información), dado que divide el continente americano en américa del sur, del norte y central. Pero notar que no todos las exportaciones que van con destino Alemania, por ejemplo, deben necesariamente desembarcar en Europa. Probablemente algunas desembarquen en Asia. Pero si, es de esperar, que la gran mayoría desembarque en Alemania. Esto será explorado mas adelante).
Países de destino
Notar que en el siguiente elemento visual, el área del cuadrado es proporcional al número de exportaciones destinadas a cada país respecto del total de exportaciones para los 2.5mm de registros existentes (no es muestral a diferencia de los otros elementos visuales explorados previamente).
datos[, .(.N), by = .(ingles, zona_geografica_puerto_desembarque)] %>%
ggplot(aes(area = N, fill= zona_geografica_puerto_desembarque, label= ingles,
subgroup=zona_geografica_puerto_desembarque)) +
geom_treemap() +
guides(fill = guide_legend(title = "Zona desembarque", title.position = "top")) +
geom_treemap_text(fontface = "italic", colour = "white", place = "centre") +
scale_fill_tableau()La importancia de USA, de las principales economías de América del sur, en particular Perú, Japón y China son bastante importantes, México, Bolivia y Paraguay, y en menor medida las principales potencias europeas.
Atípico caso de Paraguay y Bolivia
## zona_geografica_puerto_desembarque N
## 1: América Central 181536
## 2: América del Sur 2715
## 3: América del Norte 2115
## 4: Asia 84
## 5: Europa 66
## 6: Africa 27
## 7: Oceania 13
La gran mayoría de las mercaderías con destino Paraguay o Bolivia desembarcan en América Central
datos[ingles == "Bolivia" | ingles == "Paraguay", .(.N),
.(zona_geografica_puerto_desembarque, tipo_puerto_embarque)][order(-N)]## zona_geografica_puerto_desembarque tipo_puerto_embarque N
## 1: América Central Avanzada fronteriza 155170
## 2: América Central Aeropuerto 24379
## 3: América Central Puerto marítimo 1987
## 4: América del Sur Avanzada fronteriza 1941
## 5: América del Norte Avanzada fronteriza 1085
## 6: América del Norte Aeropuerto 1029
## 7: América del Sur Aeropuerto 530
## 8: América del Sur Puerto marítimo 244
## 9: Asia Avanzada fronteriza 81
## 10: Europa Aeropuerto 35
## 11: Europa Avanzada fronteriza 30
## 12: Africa Avanzada fronteriza 27
## 13: Oceania Avanzada fronteriza 13
## 14: Asia Aeropuerto 3
## 15: Europa Puerto marítimo 1
## 16: América del Norte Puerto marítimo 1
Es entendible que tanto Bolivia como Paraguay, al no tener acceso al mar, las zona de desembarque podría ser distinta, pero desmenuzando los datos por tipo de puerto de embarque, resutla que la gran mayoría de las exportaciones ocurrieron por avanzada fronteriza y zonas de desembarque de América Central, lo que parece cuanto menos llamativo.
Es avanzada fronteriza siempre por vía terrestre?
## glosa_viatransporte tipo_puerto_embarque N
## 1: CARRETERO / TERRESTRE Avanzada fronteriza 158295
## 2: FERROVIARIO Avanzada fronteriza 51
## 3: OTRA Avanzada fronteriza 1
Avanzada fronteriza es practicamente siempre carretero terrestre. Pero pareciera demasiado poco eficiente o un poco dificil de digerir que las exportaciones vayan a las Antillas Holandesas para que luego viajen en sentido contrario a Bolivia o Paraguay.
Mirando en mayor detalle el caso…
bol_par <- (datos[ingles == "Bolivia" | ingles == "Paraguay", .(.N), .(zona_geografica_puerto_desembarque,
tipo_puerto_embarque,
pais_puerto_desembarque)])[order(-N)] %>% filter(N>30)
bol_par %>%
ggplot() +
geom_histogram(aes(fct_reorder(pais_puerto_desembarque, N), N, fill = zona_geografica_puerto_desembarque),
stat = 'identity') +
theme_bw() +
coord_flip() +
theme(axis.title.y = element_blank()) +
scale_fill_tableau() +
guides(fill = guide_legend(title = "")) +
facet_wrap(~tipo_puerto_embarque) +
labs(title = "Por cual vía de transporte y a que país desembarca\nlas mercaderías que tienen como destino Bolivia y Paraguay") +
theme(plot.title = element_text(face = "bold"))## [1] 0.9735899
Un 97.2% de las exportaciones que van dirigidas a Bolivia o Paraguay, que son mas de 180000, pasaron por las Antillas Holandesas. De las cuales:
sum(bol_par[pais_puerto_desembarque == "Antillas Holandesas" &
tipo_puerto_embarque == "Avanzada fronteriza",
N])/sum(bol_par[pais_puerto_desembarque == "Antillas Holandesas", N])## [1] 0.8548804
Un 85% llega por vía terrestre.
Es decir, las exportaciones llegan a las antillas holandesas vía avanzada fronteriza, y luego
datos[nombre_puerto_desembarque == "OTROS PUERTOS DE AMÉRICA NO ESPECIFICADOS",
.(.N),.(tipo_puerto_embarque)]## tipo_puerto_embarque N
## 1: Aeropuerto 42767
## 2: Puerto marítimo 153642
## 3: Avanzada fronteriza 159141
A las antillas Holandesas llegan exportaciones por las 3 vías. Y a pesar de que las distancia a las Antillas Holandesas no es poca, hay muchísimo que llega por via terrestre, y mayor curiosidad, a destinos mucho mas cercanos, como el caso de Paraguay o Bolivia.
Zoom puesto en las Antillas Holandesas…
Lo único realmente atípico es el caso Bolivia y Paraguay, dado que el resto de los países de destino de mercadería que llega a las Antillas Holandesas son países centroamericanos, con unos casos marginales de Brasil, Uruguay , Argentina y Perú.
## año mes item_sa_operacionexpo valor
## Min. :2010 Min. : 1.000 Min. : 30000 Min. : -233600
## 1st Qu.:2012 1st Qu.: 3.000 1st Qu.: 8119090 1st Qu.: 3129
## Median :2015 Median : 6.000 Median :22042168 Median : 17346
## Mean :2015 Mean : 6.169 Mean :34640395 Mean : 256996
## 3rd Qu.:2017 3rd Qu.: 9.000 3rd Qu.:48171000 3rd Qu.: 66484
## Max. :2020 Max. :12.000 Max. :99999999 Max. :1208175003
##
## peso long lat pib
## Min. : 0 Min. :-177.39 Min. :-85.25 Min. : 210.8
## 1st Qu.: 450 1st Qu.: 0.92 1st Qu.: 4.59 1st Qu.: 6453.6
## Median : 5751 Median : 32.65 Median : 14.90 Median : 11079.7
## Mean : 225973 Mean : 45.08 Mean : 11.86 Mean : 25369.8
## 3rd Qu.: 26844 3rd Qu.: 105.82 3rd Qu.: 26.87 3rd Qu.: 48919.8
## Max. :1084713000 Max. : 178.24 Max. : 77.96 Max. :110701.9
## NA's :56563 NA's :56563 NA's :89109
## lluvia lval lpes
## Min. : 56 Min. :-Inf Min. : 0.000
## 1st Qu.: 715 1st Qu.: 8 1st Qu.: 6.000
## Median :1130 Median : 9 Median : 8.000
## Mean :1286 Mean :-Inf Mean : 7.619
## 3rd Qu.:1738 3rd Qu.: 11 3rd Qu.:10.000
## Max. :3240 Max. : 20 Max. :20.000
## NA's :112049 NA's :983
La variable fob_us_inv tiene valores negativos.
## [1] 983
Mirando en mayor detalle estas 2 variables que tienen prácticamente la misma media y la distribución por cuartil sigue evolución similar. Además el valor de la exportación es de particular interés.
select_cols = c("valor", "peso")
datos[ , ..select_cols] %>% apply(., 2, function(x) quantile(x, probs = seq(0,1,0.1)))## valor peso
## 0% -233600.000 0.370
## 10% 602.196 33.200
## 20% 2048.200 217.000
## 30% 4600.000 834.000
## 40% 9312.410 2250.408
## 50% 17345.630 5751.200
## 60% 29280.000 13936.000
## 70% 49938.604 22908.000
## 80% 91634.038 43680.000
## 90% 216650.980 109785.832
## 100% 1208175003.000 1084713000.000
Como se puede ver, tanto el peso como el valor siguen aproximadamente una distribución log-normal. Esto se puede explorar en mayor profundidad, testeando a cual distribución se asemeja de mejor manera mediante el método de máxima verosimilitud, que brinda el paquete univariateML. Solo hay que hacer una pequeña corrección, dado que el paquete requiere que los valores mínimos del x sean mayores a 0, y la data tiene 983 valores negativos, y 1317 valores que estan entre 0 y 1, cuyo log es igual o menor a 0 (el logaritmo de 1 es 0). Por tanto, se hace un pequeño ajuste que genera ruido, que ensucia los datos, que es botar los datos negativos y sumar un numero marginalmente mayor a 1 a todos el resto de los valores, de modo tal que se intenta preservar la integridad de los datos originales a la vez que ajustarse a los requiermientos del método (el método no puede trabajar con datos que tengan probabilidad 0 o menor). Potencialmente, un mejor ajuste tendría que ser gradual, de modo que preserve la forma de la distribución. Con ello, en lugar de sumar una constante a todos los valores, se podría sumar un cantidad que varíe a lo largo de la distribución, (posiblemente) incrementando su magnitud para los valores mas altos. Esto escapa, por ahora, el enfoque del estudio.
datos. <- datos[!datos$valor < 0, ]
datos. <- na.omit(datos.)
datos.$valor <- 1.000001+ datos.$valor
comp_aic <- AIC(
mlbetapr(log(datos.$valor)),
mlexp(log(datos.$valor)),
mlinvgamma(log(datos.$valor)),
mlgamma(log(datos.$valor)),
mlnorm(log(datos.$valor)),
mllnorm(log(datos.$valor)),
mlrayleigh(log(datos.$valor)),
mlinvgauss(log(datos.$valor))
)
comp_aic## df AIC
## mlbetapr(log(datos.$valor)) 2 12893393
## mlexp(log(datos.$valor)) 1 16520275
## mlinvgamma(log(datos.$valor)) 2 26998148
## mlgamma(log(datos.$valor)) 2 11939519
## mlnorm(log(datos.$valor)) 2 11570286
## mllnorm(log(datos.$valor)) 2 12496493
## mlrayleigh(log(datos.$valor)) 1 13497824
## mlinvgauss(log(datos.$valor)) 2 32708019
Entre las distribuciones BetaPrime, exponencial, inversa de la gamma, gamma, normal, log-normal, reyleigh, y la inversa de la gaussiana, para los valores logarítmicos del valor de exportación (notar que log-normal de los valores logaritmicos es el logaritmo del logaritmo del valor), la distribución mas apropiada es la normal. Es decir, para los valores en dolar de las exportaciones, la mejor distribución es la log-normal.
datos. <- datos[!datos$peso < 0, ]
datos. <- na.omit(datos.)
datos.$peso <- 1.000001+ datos.$peso
comp_aic <- AIC(
mlbetapr(log(datos.$peso)),
mlexp(log(datos.$peso)),
mlinvgamma(log(datos.$peso)),
mlgamma(log(datos.$peso)),
mlnorm(log(datos.$peso)),
mllnorm(log(datos.$peso)),
mlrayleigh(log(datos.$peso)),
mlinvgauss(log(datos.$peso))
)
comp_aic## df AIC
## mlbetapr(log(datos.$peso)) 2 15108661
## mlexp(log(datos.$peso)) 1 15704879
## mlinvgamma(log(datos.$peso)) 2 15737715
## mlgamma(log(datos.$peso)) 2 13675179
## mlnorm(log(datos.$peso)) 2 12985273
## mllnorm(log(datos.$peso)) 2 14407871
## mlrayleigh(log(datos.$peso)) 1 13489794
## mlinvgauss(log(datos.$peso)) 2 14926802
Idem para la variable peso.
El diagnóstico más formal de semejanza con la distribución normal es mediante el test de shapiro, pero tiene el inconveniente, al igual que el resto de test de normalidad (hay argumentos muy interesantes al respecto en el foro estadístico Cross Validated) de siempre rechazar normalidad cuando las muestras toman valores importantes. Para testear esto, trabaje con una segunda muestra de 5000 y otra de 50 datos, graficaré un qqplot comparando lo teórico con lo muestral para cada muestra y un test de normalidad de shapiro.
##
## Shapiro-Wilk normality test
##
## data: log(muestra2$valor)
## W = 0.99631, p-value = 0.0000000007816
-Se rechaza abruptamente la hipótesis de que la muestra provenda de una distribución normal.
##
## Shapiro-Wilk normality test
##
## data: log(muestra2$valor)
## W = 0.97224, p-value = 0.2848
Se recomienda entonces mejor determinar si la distribución asemeja en mejor medida a la distribución normal mediante elementos visuales, principalmente el qqplot. (pero también se puede utilizar la función de distribución acumulada o comparar visualmente la distribución del valor con una normal que tenga los mismos parámetros de media y desviación estandar (este trabajo lo realicé pero agregando distintas variables, que agregaré en otro documento), o incluso comparar incorporando el elemento histograma).
Pib, lluvia e Item de operación
El caracter multimodal de las distribuciones se explica, en parte, por que se trata de data de países, que presentan claramente ausencia de aleteoridad en su número de ocurrencias. Es posible que cada país tenga una participación relativamente estable en el tiempo -esa es la hipótesis-, o que tengan bastante aleteoridad, pero en rango, de tal manera que países muy activos como USA o Perú, varíen el número de exportaciones pero se mantenga siempre en rangos altos, y otros lo contrario.
Para el calculo de las modas de las distribuciones, me apoyaré de la función creada acá, que explora la detección de múltiples modas en funciones de estas características.
modes <- function(d){
i <- which(diff(sign(diff(d$y))) < 0) + 1
data.frame(x = d$x[i], y = d$y[i])
}
ggplot(data.frame(x = densidad_pib$x, y = densidad_pib$y * densidad_pib$n), aes(x, y)) +
geom_density(stat = 'identity', fill = 'tomato2', alpha = .3) +
geom_point(data = modes(densidad_pib), aes(y = y * densidad_pib$n)) +
theme_minimal() +
labs(title = "Modas de la densidad del PIB de los países",y = "", x="")## x y
## 1 7247.232 0.00005830030984
## 2 26865.832 0.00000465570714
## 3 33564.866 0.00000601728550
## 4 54618.973 0.00003982140339
## 5 63710.519 0.00000199895992
## 6 78783.346 0.00000052181633
## 7 92181.414 0.00000059519006
## 8 110603.757 0.00000003849748
## nombre_pais PIB N
## 1: Estados Unidos 54659.200 15381
## 2: Perú 6453.561 10566
## 3: Bolivia 2559.511 5259
## 4: Argentina 10043.510 5157
## 5: Brasil 11079.710 4762
## 6: Colombia 7691.746 4487
## 7: China 7752.560 4112
## 8: México 10403.540 3660
## 9: Ecuador 5185.091 3240
## 10: Japón 48919.800 3165
## 11: Canadá 51393.220 3031
## 12: Países Bajos 55020.980 2749
## 13: Reino Unido 43343.270 2688
## 14: Uruguay 14617.460 2473
## 15: España 32898.330 2135
Los otros 2 pequeños saltos de la distribución (entre los 2 grandes), seguramente se trata de España y Corea del Sur (más otros países de menor frecuencia de compra con similar poder adquisitivo). “Segundo Mundo o países quasi-desarrollados”.
Es una función mas dificil de caracterizar, pero uno podría evaluar los 3 principales máximos/modas.
## x y
## 1 78.98767 0.00006278128
## 2 703.78453 0.00169362514
## 3 1167.11703 0.00060069375
## 4 1735.75237 0.00108484485
## 5 1995.49938 0.00012532514
## 6 2276.30695 0.00022791999
## 7 2922.16438 0.00020194616
## 8 3238.07290 0.00026537055
3 grandes grupos: 700mm, 1166mm y 1700mm.
## nombre_pais lluvia N
## 1: Estados Unidos 715 15381
## 2: Perú 1738 10566
## 3: Bolivia 1146 5259
## 4: Argentina 591 5157
## 5: Brasil 1761 4762
## 6: Colombia 3240 4487
## 7: China 645 4112
## 8: México 758 3660
## 9: Ecuador 2274 3240
## 10: Japón 1668 3165
## 11: Canadá 537 3031
## 12: Países Bajos 778 2749
## 13: Reino Unido 1220 2688
## 14: Uruguay 1300 2473
## 15: España 636 2135
Claramente no es data tan intersante dado que en toda descomposición simplemente vemos reflejados las variables en los principales compradores. Por ejemplo, el máximo es nuevamente USA compartiendo grupo con China, Argentina, México Canada y Páises bajos, luego Reino Unido y Bolivia y finalmente países igualmente lluviosos como Perú y Brasil.
Es evidente que este tipo de análisis no es revelador de información sensible, sino el reflejo de cuanto en común para cada variable comparten los principales países. Pero, agregando variables de estacionalidad y del tipo de la mercadería que exportan (desglosando en tiempo), este tipo de información puede ayudar a generar pequeños ajustes a modelos predictivos (de estimación de demanda por tipo de producto en el tiempo por ejemplo).
La variable item es independiente de la variable país de destino. Trata de una categorización desconocida de (probablemente) la mercadería, y es interesante estudiar que tipo de información revelan las modas de la función de densidad, pero dado que esto es un análisis exploratorio inicial, y que se requiere interactuar con otras variables -de momento desconocidas- para entender mejor el caso, entonces solo dejaré expuesto el mismo procedimiento anterior.
## x y
## 1 7887528 0.000000031552584
## 2 21839172 0.000000036069970
## 3 30840232 0.000000007670695
## 4 40516373 0.000000011749142
## 5 62118918 0.000000003307146
## 6 73145217 0.000000007676519
## 7 84846596 0.000000018073435
Las densidades de las modas son significativamente menores a la de los casos anteriores, insinuando que tal vez haya mucho ruido en la variable.
En concreto, una muestra de 5000 registros será explorada para entender la magnitud de la cantidad de información, ruido o entropía que contiene (en términos de la teoría de la información):
set.seed(121)
muestra_peq <- muestra[sample(nrow(muestra), 5000),]
a1 <- muestra_peq %>% data.frame(item = .$item_sa_operacionexpo, indice = seq_along(.$item_sa_operacionexpo)) %>% ggplot() + geom_point(aes(item, indice)) + theme_tufte()
a2 <- muestra_peq %>% data.frame(item = log(.$item_sa_operacionexpo), indice = seq_along(.$item_sa_operacionexpo)) %>% ggplot() + geom_point(aes(item, indice)) + theme_tufte()
a3 <- muestra_peq %>% data.frame(item = sqrt(.$item_sa_operacionexpo), indice = seq_along(.$item_sa_operacionexpo)) %>% ggplot() + geom_point(aes(item, indice)) + theme_tufte()
a. <- ggarrange(a1, a3, a2)
title_. <- expression(atop(bold("Distribución exacta variable Item"),
scriptstyle("Escala normal, logarítmica y cuadrática\nTamaño Muestral = 5000")))
annotate_figure(a., top = text_grob(title_.))## 0% 10% 20% 30% 40% 50% 60% 70%
## 198900 8055010 8093010 17041066 22042141 22042168 33049910 44091022
## 80% 90% 100%
## 70134900 84439990 99999999
Y la función de distribución acumulada para
ggplot(muestra_peq,aes(item_sa_operacionexpo)) +
stat_ecdf(geom="point") +
labs(title = "") + theme_minimal()## nombre_pais item N
## 1: Estados Unidos 30415745 15381
## 2: Perú 54257667 10566
## 3: Bolivia 56325222 5259
## 4: Argentina 53087848 5157
## 5: Brasil 31159153 4762
## 6: Colombia 38837017 4487
## 7: China 22476644 4112
## 8: México 36536767 3660
## 9: Ecuador 40699217 3240
## 10: Japón 18936716 3165
## 11: Canadá 26579561 3031
## 12: Países Bajos 22114559 2749
## 13: Reino Unido 20304459 2688
## 14: Uruguay 50456330 2473
## 15: España 23819947 2135
## 16: Paraguay 45179479 2090
## 17: Alemania 28853765 1845
## 18: Corea del Sur 21709872 1786
## 19: Costa Rica 33694758 1753
## 20: Panamá 32498305 1339
Puede que esta variable contenga un ordenamiento lógico, dado que países con mucho movimiento como USA, Perú o Japón tienen un valor de item promedio bastante disímil. Y algunos países tienen mucha similitud como Argentina-Perú-Bolivia o China-Canada-Países Bajos-reino Unido. Es interesante mirar si los promedios por continente entregan alguna pista al respecto.
## zona_geografica_puerto_desembarque item N
## 1: América del Sur 46375950 30678
## 2: América del Norte 30019180 18972
## 3: América Central 40919329 18772
## 4: Europa 24235401 15934
## 5: Asia 22034813 13424
## 6: Oceania 35380778 1334
## 7: Africa 31988061 886
Efectivamente parece haber un agrupamiento por zona geográfica, tal como se vio por país. El orden pareciera ser: - Asia, Europa - Oceania, Africa, America del Norte, Oceanía - America Central, América del Sur
## zona_geografica_puerto_desembarque estac item N
## 1: Africa inv 35261129 200
## 2: Africa pri 27336775 254
## 3: Africa ver 32361907 251
## 4: Africa oto 34380204 181
## 5: América Central ver 41382662 4496
## 6: América Central inv 39913678 4796
## 7: América Central oto 40316203 5009
## 8: América Central pri 42207858 4471
## 9: América del Norte pri 28229553 5598
## 10: América del Norte ver 31179342 4665
## 11: América del Norte inv 28470901 5107
## 12: América del Norte oto 33493148 3602
## 13: América del Sur pri 46357758 7331
## 14: América del Sur ver 46136891 7469
## 15: América del Sur inv 47169415 7792
## 16: América del Sur oto 45848647 8086
## 17: Asia ver 22700958 3201
## 18: Asia pri 20554634 4068
## 19: Asia oto 23903624 2663
## 20: Asia inv 21723358 3492
## 21: Europa ver 24424568 4166
## 22: Europa pri 20968949 5172
## 23: Europa oto 28051656 2892
## 24: Europa inv 25604032 3704
## 25: Oceania oto 35760520 392
## 26: Oceania pri 34178075 306
## 27: Oceania ver 34500676 274
## 28: Oceania inv 36652371 362
## zona_geografica_puerto_desembarque estac item N
Y la tendencia se sostiene a lo largo de las estaciones. Realizando este reordenamiento uno puede aprovechar de notar cuanto efecto tiene el otoño estadounidense en el número de exportaciones que realizan.
Para tener una noción visual al respecto:
datos[, .(item = mean(item_sa_operacionexpo), .N),
.(zona_geografica_puerto_desembarque, estac)][order(zona_geografica_puerto_desembarque)] %>%
ggplot(aes(x=fct_reorder(estac, item), y = item, fill = item)) +
geom_bar(colour="black", stat="identity",
position=position_dodge(), size = .3) +
theme_bw() +
guides(fill = 'none') +
labs(title = "La 'categorización' de los ítem se sostiene a lo largo de las estaciones",
subtitle = "Pero de modo dispar entre las zonas geográficas",
x = "", y = "Item promedio") +
facet_wrap(~zona_geografica_puerto_desembarque, ncol = 4) +
scale_x_discrete(label = function(x) stringr::str_sub(x, -27)) +
theme(plot.title = element_text(size = 14))Aprovechando que el análisis de la unidad anterior levantó cierta estacionalidad (en particular la reducción de USA en otoño), se expande un poco mas al respecto.
Notar como en Otoño se encogen las “compras” de América del Norte, Europa y Asia (y África en menor medida), y como estas se expanden en Primavera. Curiosamente en Verano es cuando mas se encogen las “compras” realizadas por países sudamericanos. (recordar que la variable es la estación del país de destino). Esta relación se explorará con mayor detenimiento más adelante.
Notar que la excepcionalidad de otoño en América del Norte también esta presente en el gráfico anterior a este, donde se puede ver que este subcontinente aumenta el item promedio, más cercano a los de América del Sur y Central. ¿Tendrá el item cierta relación con el valor?. Notar que Asia y Europa (también en el gráfico anterior a este) tienen valores promedios de Item mas bajos. Estas relaciones se explorarán mas adelante.
También queda patente que en los otoños de Africa, Asia, Europa y América del Norte, los promedios de item son los máximos y los de primavera son los mínimos.
Obviamente el año 2020 esta incompleto, es el único?…
## # A tibble: 11 x 2
## año n
## <int> <int>
## 1 2010 9
## 2 2011 12
## 3 2012 12
## 4 2013 12
## 5 2014 12
## 6 2015 12
## 7 2016 12
## 8 2017 12
## 9 2018 12
## 10 2019 12
## 11 2020 4
## # A tibble: 12 x 2
## mes n
## <int> <int>
## 1 1 11
## 2 2 11
## 3 3 11
## 4 4 11
## 5 5 10
## 6 6 10
## 7 7 10
## 8 8 10
## 9 9 10
## 10 10 9
## 11 11 9
## 12 12 9
Las implicancias es que los meses estan desbalanceados, con 4 meses disponibles para los 11 años, 5 para 10, 3 para solo 9.
De momento, quitando 2010 y 2020.
Claramente hay factor de temporalidad reflejado en la variabilidad de los meses, lo cual tendrá que estudiarse en mayor detalles más adelante.
muestra[año != 2010 & año != 2020,] %>%
ggplot() +
geom_line(aes(as.yearmon(fecha), group = 1), stat = "count", color = "lightblue", alpha = 2, size = 2) +
theme_tufte() +
theme(axis.text.x = element_text(angle = 90, size = 10)) +
scale_y_continuous(limits = c(500, 1100)) +
labs(x = "", title = "Número de exportaciones mensuales",
subtitle = "Ocurridas entre 2011 y 2019")La variable año-mes muestra gran variabilidad en el tiempo, empujado fundamentalmente por la variable mes, que es la mas ruidosa.
Adicionalmente, agregue a los datos la variable estación, ya que luego de obtener los datos de longitud y latitud, se puede establecer de modo simple que aquellos que estan sobre la línea del ecuador tienen cierta estación del año distinta a las que estan por bajo la línea del ecuador, aunque hay ciertas diferencias intra-linea ecuatorial, como las diferencias entre China e India, u otros donde el inicio de las estaciones es un poco mas tardío en algunos que en otros.