EDA dataset R

El dataset fue conseguido del UCI en el link https://archive.ics.uci.edu/dataset/585/stock+keeping+units, esta abierto para descarga de cualquiera

df <- read_excel("stock+keeping+units/sku_data.xlsx") %>%
  janitor::clean_names()  # Limpia nombres: sin espacios, minúsculas y con "_"

Descripcion de las variables

1. unit_price - unit price in euro

2. expire_date - shelf-life

3. total_outbound - number of pallets sold from 2017-02-06 to 2018-02-13

4. outbound_number - how many times a product was ordered from 2017-02-06 to 2018-02-13

5. pal_grossweight - how much a fully-loaded pallet weights (kg)

6. pal_height - height of a fully-loaded pallet (cm)

7. units_per_pal - units per pallet

head(df)

## # A tibble: 6 × 8
##      id unitprice expire_date outbound_number total_outbound pal_grossweight
##   <dbl>     <dbl>       <dbl>           <dbl>          <dbl>           <dbl>
## 1     1     0.058         547               9           2441            106.
## 2     2     0.954         547               0              0            208.
## 3     3     2.38          547              12             23            166.
## 4     4     5.1           547               0              0            221.
## 5     5     0             547               0              0              0 
## 6     6     1.11          547               1              1            208.
## # ℹ 2 more variables: pal_height <dbl>, units_per_pal <dbl>

str(df)

## tibble [2,279 × 8] (S3: tbl_df/tbl/data.frame)
##  $ id             : num [1:2279] 1 2 3 4 5 6 7 8 9 10 ...
##  $ unitprice      : num [1:2279] 0.058 0.954 2.385 5.1 0 ...
##  $ expire_date    : num [1:2279] 547 547 547 547 547 547 547 547 547 365 ...
##  $ outbound_number: num [1:2279] 9 0 12 0 0 1 2 0 0 0 ...
##  $ total_outbound : num [1:2279] 2441 0 23 0 0 ...
##  $ pal_grossweight: num [1:2279] 106 208 166 221 0 ...
##  $ pal_height     : num [1:2279] 1.56 1 1.02 1.05 0 1 1.02 0 0 1.44 ...
##  $ units_per_pal  : num [1:2279] 1920 384 108 72 0 384 108 0 0 36 ...

glimpse(df)

## Rows: 2,279
## Columns: 8
## $ id              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
## $ unitprice       <dbl> 0.0580, 0.9540, 2.3850, 5.1000, 0.0000, 1.1100, 2.5150…
## $ expire_date     <dbl> 547, 547, 547, 547, 547, 547, 547, 547, 547, 365, 547,…
## $ outbound_number <dbl> 9, 0, 12, 0, 0, 1, 2, 0, 0, 0, 0, 6, 8, 0, 0, 1, 0, 3,…
## $ total_outbound  <dbl> 2441, 0, 23, 0, 0, 1, 4, 0, 0, 0, 0, 14, 12, 0, 0, 2, …
## $ pal_grossweight <dbl> 105.60, 207.68, 165.78, 221.04, 0.00, 207.68, 165.78, …
## $ pal_height      <dbl> 1.56, 1.00, 1.02, 1.05, 0.00, 1.00, 1.02, 0.00, 0.00, …
## $ units_per_pal   <dbl> 1920, 384, 108, 72, 0, 384, 108, 0, 0, 36, 72, 384, 10…

dim(df) # Filas y columnas

## [1] 2279    8

summary(df)

##        id           unitprice        expire_date    outbound_number
##  Min.   :   1.0   Min.   :  0.000   Min.   :  0.0   Min.   :   0   
##  1st Qu.: 570.5   1st Qu.:  0.000   1st Qu.:365.0   1st Qu.:   0   
##  Median :1140.0   Median :  1.294   Median :547.0   Median :   1   
##  Mean   :1140.0   Mean   :  4.269   Mean   :410.4   Mean   : 236   
##  3rd Qu.:1709.5   3rd Qu.:  4.545   3rd Qu.:547.0   3rd Qu.:  45   
##  Max.   :2279.0   Max.   :518.592   Max.   :734.0   Max.   :6325   
##  total_outbound    pal_grossweight   pal_height     units_per_pal     
##  Min.   :    0.0   Min.   :  0.0   Min.   :0.0000   Min.   :     0.0  
##  1st Qu.:    0.0   1st Qu.: 60.0   1st Qu.:0.0000   1st Qu.:    32.0  
##  Median :    3.0   Median :167.7   Median :0.8400   Median :   108.0  
##  Mean   :  731.7   Mean   :192.9   Mean   :0.6728   Mean   :   755.6  
##  3rd Qu.:  419.5   3rd Qu.:277.6   3rd Qu.:1.0200   3rd Qu.:   384.0  
##  Max.   :26411.0   Max.   :907.2   Max.   :2.1600   Max.   :200000.0

skim(df) # Más detallado

Data summary

Name	df
Number of rows	2279
Number of columns	8
_______________________
Column type frequency:
numeric	8
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
id	0	1	1140.00	658.03	1	570.5	1140.00	1709.50	2279.00	▇▇▇▇▇
unitprice	0	1	4.27	14.45	0	0.0	1.29	4.54	518.59	▇▁▁▁▁
expire_date	0	1	410.37	240.88	0	365.0	547.00	547.00	734.00	▃▁▂▇▁
outbound_number	0	1	235.98	700.23	0	0.0	1.00	45.00	6325.00	▇▁▁▁▁
total_outbound	0	1	731.70	2146.03	0	0.0	3.00	419.50	26411.00	▇▁▁▁▁
pal_grossweight	0	1	192.94	164.62	0	60.0	167.68	277.56	907.20	▇▆▁▁▁
pal_height	0	1	0.67	0.55	0	0.0	0.84	1.02	2.16	▇▂▇▂▁
units_per_pal	0	1	755.56	6278.44	0	32.0	108.00	384.00	200000.00	▇▁▁▁▁

# Duplicados
sum(duplicated(df))

## [1] 0

colSums(is.na(df))

##              id       unitprice     expire_date outbound_number  total_outbound 
##               0               0               0               0               0 
## pal_grossweight      pal_height   units_per_pal 
##               0               0               0

plot_missing(df) # Gráfico de valores perdidos (de DataExplorer)

Inicialmente parece que el dataset no tiene faltantes

colSums(df == 0, na.rm = TRUE)

##              id       unitprice     expire_date outbound_number  total_outbound 
##               0             710             529             979             979 
## pal_grossweight      pal_height   units_per_pal 
##             371             787             295

colMeans(df == 0, na.rm = TRUE) # proporción

##              id       unitprice     expire_date outbound_number  total_outbound 
##       0.0000000       0.3115401       0.2321194       0.4295744       0.4295744 
## pal_grossweight      pal_height   units_per_pal 
##       0.1627907       0.3453269       0.1294427

Pero si observamos la proporcion de 0’s que tiene el dataset, estas son inusualmente altas para todas las columnas importantes. Ademas por el contexto de las variables, ninguna de estas tienen sentido que esten igualadas a 0, por lo que lo mas seguro es que el sistema que haya coleccionado los datos haya sido codificado para introducir datos faltantes como “0”, por lo que tendremos que tratar los 0’s como faltantes.

table(df$expire_date)

## 
##    0  365  456  547  730  734 
##  529  305    8 1250  186    1

Otra cosa a notar es la variable expire_date, la cual a pesar de ser tratada como numerica en el dataset, tiene una cantidad limitada de variables unicas, lo cual hace que se comporte mas como una variable categorica ordinal. Esto complicara un poquito la imputacion de datos faltantes, por lo cual veremos mas adelante como lidiaremos con esto

# Transformamos 0's a NA's
df <- df %>% mutate(across(everything(), ~ na_if(., 0)))

summary(df)

##        id           unitprice         expire_date    outbound_number 
##  Min.   :   1.0   Min.   :  0.0033   Min.   :365.0   Min.   :   1.0  
##  1st Qu.: 570.5   1st Qu.:  1.0600   1st Qu.:547.0   1st Qu.:   3.0  
##  Median :1140.0   Median :  2.8000   Median :547.0   Median :  17.0  
##  Mean   :1140.0   Mean   :  6.2014   Mean   :534.4   Mean   : 413.7  
##  3rd Qu.:1709.5   3rd Qu.:  6.5460   3rd Qu.:547.0   3rd Qu.: 374.0  
##  Max.   :2279.0   Max.   :518.5920   Max.   :734.0   Max.   :6325.0  
##                   NA's   :710        NA's   :529     NA's   :979     
##  total_outbound  pal_grossweight   pal_height    units_per_pal     
##  Min.   :    1   Min.   :  0.8   Min.   :0.040   Min.   :     1.0  
##  1st Qu.:    9   1st Qu.:128.0   1st Qu.:0.900   1st Qu.:    42.0  
##  Median :  269   Median :207.7   Median :1.020   Median :   108.0  
##  Mean   : 1283   Mean   :230.5   Mean   :1.028   Mean   :   867.9  
##  3rd Qu.: 1372   3rd Qu.:293.8   3rd Qu.:1.050   3rd Qu.:   384.0  
##  Max.   :26411   Max.   :907.2   Max.   :2.160   Max.   :200000.0  
##  NA's   :979     NA's   :371     NA's   :787     NA's   :295

skim(df)

Data summary

Name	df
Number of rows	2279
Number of columns	8
_______________________
Column type frequency:
numeric	8
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
id	0	1.00	1140.00	658.03	1.00	570.50	1140.00	1709.50	2279.00	▇▇▇▇▇
unitprice	710	0.69	6.20	17.07	0.00	1.06	2.80	6.55	518.59	▇▁▁▁▁
expire_date	529	0.77	534.42	96.11	365.00	547.00	547.00	547.00	734.00	▂▁▇▁▁
outbound_number	979	0.57	413.68	886.73	1.00	3.00	17.00	374.00	6325.00	▇▁▁▁▁
total_outbound	979	0.57	1282.73	2714.59	1.00	9.00	269.00	1372.00	26411.00	▇▁▁▁▁
pal_grossweight	371	0.84	230.46	154.01	0.80	128.00	207.68	293.76	907.20	▇▇▁▁▁
pal_height	787	0.65	1.03	0.32	0.04	0.90	1.02	1.05	2.16	▁▂▇▂▁
units_per_pal	295	0.87	867.91	6722.00	1.00	42.00	108.00	384.00	200000.00	▇▁▁▁▁

plot_missing(df)

Ahora si podemos notar que este dataset tiene una cantidad considerable de datos faltantes que tendremos que afrontar usando tecnicas avanzadas de prediccion

# Histogramas 
df %>% select(where(is.numeric), -id) %>% 
  gather(variable, valor) %>%
  ggplot(aes(valor)) +
  geom_histogram(bins = 30, fill = "skyblue") +
  facet_wrap(~variable, scales = "free")

# Boxplot para ver outliers
df %>%
  select(where(is.numeric), -id) %>%  
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = "", y = valor)) +
  geom_boxplot(fill = "red2") +
  facet_wrap(~variable, scales = "free_y") +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())

# Hacemos pruebas de normalidad
df %>%
  select(where(is.numeric)) %>%
  summarise(across(everything(), 
                   ~ shapiro.test(.x)$p.value))

## # A tibble: 1 × 8
##         id unitprice expire_date outbound_number total_outbound pal_grossweight
##      <dbl>     <dbl>       <dbl>           <dbl>          <dbl>           <dbl>
## 1 5.59e-26  7.49e-63    7.20e-48        3.34e-50       3.97e-51        2.51e-34
## # ℹ 2 more variables: pal_height <dbl>, units_per_pal <dbl>

# Ninguna variable es normal 

Tratamiento de atipicos

Vamos a inicialmente tratar los atipicos, para luego tratarlos con los faltantes que ya tenemos.

Por la mayor parte pareceria que varias columnas tienen varios datos faltantes que deberiamos de eliminar, pero como estos son datos reales no es tan necesario eliminarlos a todos, de todas formas hay columnas con unos datos un poquito exagerados que si creeremos necesarios corregir. Entonces esto es lo que haremos

• id: Nada
• expire_date: Por lo ya descrito, no tendria sentido tratar atipicos en esta.
• total_outbound y outbound_number: Por el contexto de las variables, estas dos estan correlacionadas, por lo que podriamos perder informacion valiosa si decidimos imputar atipicos sin tener cuidado. De todas formas viendo el boxplot, los datos alejados de su mediana no parecen tanto como anomalias sino como una tendencia que tiene sentido debido a su naturaleza, entonces estas dos variables las dejaremos quietas.
• units_per_pal y unit_price: Como podemos ver en el boxplot, hay ciertas instancias con valores muy aislados de como se comportan el resto de observaciones, a pesar de que estas sean cifras reales, vemos la necesidad de corregir estas exageraciones. Para esto vamos a aplicar el test de Grubbs tipo 10, buscando analizar si el valor maximo es atipico, esto lo hacemos iterativamente hasta que el test no encuentre mas valores aislados.
• pal_height y pal_grossweight: No deberian ser muy significantes sus atipicos, pero aun asi vamos a aplicar percentiles de 1 y 99 para corregir usando capping.

# Funcion para hacer capping usando percentiles
cap_percentiles <- function(x, lower_p = 0.01, upper_p = 0.99) {
  lower <- quantile(x, lower_p, na.rm = TRUE)
  upper <- quantile(x, upper_p, na.rm = TRUE)
  
  x <- ifelse(x < lower, lower, x)
  x <- ifelse(x > upper, upper, x)
  return(x)
}

# Función para aplicar Grubbs tipo 10 iterativamente
remove_grubbs_max <- function(x, alpha = 0.01) {
  repeat {
    test <- grubbs.test(x, type = 10)
    if (test$p.value < alpha) {
      max_val <- max(x, na.rm = TRUE)
      x[which.max(x)] <- NA   # Reemplaza solo uno por NA
    } else {
      break
    }
  }
  return(x)
}

df <- df %>%
  mutate(across(c(pal_height, pal_grossweight), cap_percentiles))

df <- df %>%
  mutate(across(c(units_per_pal, unitprice), remove_grubbs_max))

# Hacemos los boxplots de nuevo para ver cambios
df %>%
  select(where(is.numeric), -id) %>%  
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = "", y = valor)) +
  geom_boxplot(fill = "red2") +
  facet_wrap(~variable, scales = "free_y") +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())

Datos faltantes

plot_missing(df) 

Ahora para los datos faltantes, las variables de units_per_pal y pal_gross_weight no tienen porcentajes muy altos por lo que imputaremos por mediana. Mientras tanto, el resto de variables si tienen porcentajes que superan el 30% por lo que aplicaremos mice con metodo pmm para imputarlos, pero por lo descrito anteriormente del expire_date preferiremos imputar esta columna especificamente con knn, para que el modelo no se invente valores no existentes. El knn sera aplicado despues del mice para tener mejores resultados

# Imputamos por mediana por units_per_pal y pal_gross_weight
df <- df %>%
  mutate(across(c(units_per_pal, pal_grossweight), ~ ifelse(is.na(.), median(., na.rm = TRUE), .)))


cols_mice <- c("unitprice", "pal_height", "total_outbound", "outbound_number") 
col_knn   <- "expire_date"

# Guardamos columnas originales para hacer comparaciones mas adelante
df_original <- df %>% select(all_of(c(cols_mice, col_knn)))

# Guardar versión sin NA's para graficar
df_no_na <- lapply(df_original, function(col) col[!is.na(col)])

# MICE
mice_model <- mice(df[cols_mice], method = "pmm", m = 1, seed = 123)
df[cols_mice] <- complete(mice_model)

# KNN
df[col_knn] <- kNN(df, variable = "expire_date", k = 5, imp_var = FALSE)[, col_knn]

plot_missing(df)

Ahora validamos normalidad de nuevo

# df_original es antes de imputacion, df es despues de imputacion
for (col in c(cols_mice, col_knn)) {
  cat("\nColumna:", col, "\n")
  print(shapiro.test(na.omit(df_original[[col]])))
  print(shapiro.test(df[[col]]))
}

## 
## Columna: unitprice 
## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(df_original[[col]])
## W = 0.77867, p-value < 2.2e-16
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  df[[col]]
## W = 0.77932, p-value < 2.2e-16
## 
## 
## Columna: pal_height 
## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(df_original[[col]])
## W = 0.90865, p-value < 2.2e-16
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  df[[col]]
## W = 0.90957, p-value < 2.2e-16
## 
## 
## Columna: total_outbound 
## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(df_original[[col]])
## W = 0.50049, p-value < 2.2e-16
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  df[[col]]
## W = 0.495, p-value < 2.2e-16
## 
## 
## Columna: outbound_number 
## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(df_original[[col]])
## W = 0.52808, p-value < 2.2e-16
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  df[[col]]
## W = 0.52164, p-value < 2.2e-16
## 
## 
## Columna: expire_date 
## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(df_original[[col]])
## W = 0.70859, p-value < 2.2e-16
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  df[[col]]
## W = 0.66696, p-value < 2.2e-16

Todas las columnas siguen siendo no normales, por lo que haremos pruebas no parametricas para validar la imputacion. Tambien generaremos graficas para comparar

# Graficar distribuciones
for (col in c(cols_mice, col_knn)) {
  temp <- data.frame(
    Valor = c(na.omit(df_original[[col]]), df[[col]]),
    Estado = rep(c("Antes", "Después"), c(length(na.omit(df_original[[col]])), nrow(df)))
  )
  
  p <- ggplot(temp, aes(x = Valor, fill = Estado)) +
    geom_density(alpha = 0.5) +
    labs(title = paste("Distribución de", col, "antes y después de imputar")) +
    theme_minimal()
  
  print(p)
}

# Prueba no parametrica de wilcoxon test
for (col in c(cols_mice, col_knn)){
  before <- df_no_na[[col]]
  after <- df[[col]]
  test_np <- wilcox.test(before, after)
  cat("\nColumna:", col, "\n")
  cat("Wilcoxon p:", test_np$p.value, "\n")
  }

## 
## Columna: unitprice 
## Wilcoxon p: 0.6549202 
## 
## Columna: pal_height 
## Wilcoxon p: 0.289961 
## 
## Columna: total_outbound 
## Wilcoxon p: 0.6311904 
## 
## Columna: outbound_number 
## Wilcoxon p: 0.6863162 
## 
## Columna: expire_date 
## Wilcoxon p: 0.1875679

De aqui podria parecer que hay algunos p valores dudosos, pero son los suficientes como para afirmar que nuestra imputacion conserva las distribuciones.

# Convertimos expire_date a caracteres para analizarlo asi
df <- df %>%
  mutate(expire_date = as.character(expire_date))

Analisis univariado

# Excluimos la columna de los ids de nuestras graficas
cols_to_analyze <- df %>% select(-id)

# Iterar sobre columnas numéricas
numeric_cols <- names(cols_to_analyze)[sapply(cols_to_analyze, is.numeric)]

for (col in numeric_cols) {
  cat("\n===== Variable:", col, "=====\n")
  
  # Calcular estadísticos
  sk <- skewness(df[[col]], na.rm = TRUE)
  kt <- kurtosis(df[[col]], na.rm = TRUE)
  
  cat("Skewness:", sk, "\n")
  cat("Kurtosis:", kt, "\n")
  
  # Interpretación:
  # Skewness: mide la asimetría de la distribución
  #   0 ≈ simétrica, >0 sesgo a la derecha, <0 sesgo a la izquierda
  # Kurtosis: mide la "apuntamiento" de la distribución
  #   3 ≈ mesocúrtica (normal), >3 leptocúrtica (picos altos), <3 platicúrtica (aplanada)
  
  # Histograma + densidad
  p <- ggplot(df, aes(x = .data[[col]])) +
    geom_histogram(aes(y = ..density..), bins = 30, fill = "skyblue", alpha = 0.6) +
    geom_density(color = "red", size = 1) +
    labs(title = paste("Distribución de", col),
         x = col, y = "Densidad") +
    theme_minimal()
  print(p)
  
}

## 
## ===== Variable: unitprice =====
## Skewness: 1.506552 
## Kurtosis: 1.351257

## 
## ===== Variable: outbound_number =====
## Skewness: 3.273482 
## Kurtosis: 12.06157

## 
## ===== Variable: total_outbound =====
## Skewness: 4.272047 
## Kurtosis: 22.9543

## 
## ===== Variable: pal_grossweight =====
## Skewness: 1.449593 
## Kurtosis: 2.846267

## 
## ===== Variable: pal_height =====
## Skewness: 0.460275 
## Kurtosis: 0.6478598

## 
## ===== Variable: units_per_pal =====
## Skewness: 2.383801 
## Kurtosis: 5.135571

AU categoricas

cat_col <- "expire_date"

ggplot(df, aes(x = .data[[cat_col]])) +
  geom_bar(fill = "lightgreen", alpha = 0.7) +
  labs(title = paste("Frecuencias de", cat_col),
       x = cat_col, y = "Conteo") +
  theme_minimal()

# Grafico de torta
df %>%
  select(where(~is.factor(.) || is.character(.))) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  group_by(variable, valor) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(variable) %>%
  mutate(
    porcentaje = n / sum(n) * 100,
    etiqueta = paste0(round(porcentaje, 1), "%")
  ) %>%
  ggplot(aes(x = "", y = porcentaje, fill = valor)) +
  geom_col(width = 1, color = "white") +
  coord_polar(theta = "y") +
  geom_text(aes(label = etiqueta),
            position = position_stack(vjust = 0.5), size = 3) +
  facet_wrap(~variable) +
  theme_void()

Analisis bivariado

# Gráfico de correlaciones
ggcorr(
  df %>% select(all_of(numeric_cols)),
  method = c("complete.obs", "spearman"), # primero 'use', luego 'method'
  label = TRUE, label_round = 2
) +
  ggtitle("Correlación Spearman entre variables numéricas")