Empezamos cargando el dataset y haciendo un analisis basico

url <- "https://raw.githubusercontent.com/stedy/Machine-Learning-with-R-datasets/master/insurance.csv"
dbs <- read.csv(url)

dim(dbs) # Filas y columnas
## [1] 1338    7
head(dbs)
##   age    sex    bmi children smoker    region   charges
## 1  19 female 27.900        0    yes southwest 16884.924
## 2  18   male 33.770        1     no southeast  1725.552
## 3  28   male 33.000        3     no southeast  4449.462
## 4  33   male 22.705        0     no northwest 21984.471
## 5  32   male 28.880        0     no northwest  3866.855
## 6  31 female 25.740        0     no southeast  3756.622
str(dbs)
## 'data.frame':    1338 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr  "female" "male" "male" "male" ...
##  $ bmi     : num  27.9 33.8 33 22.7 28.9 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr  "yes" "no" "no" "no" ...
##  $ region  : chr  "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num  16885 1726 4449 21984 3867 ...
summary(dbs)
##       age            sex                 bmi           children    
##  Min.   :18.00   Length:1338        Min.   :15.96   Min.   :0.000  
##  1st Qu.:27.00   Class :character   1st Qu.:26.30   1st Qu.:0.000  
##  Median :39.00   Mode  :character   Median :30.40   Median :1.000  
##  Mean   :39.21                      Mean   :30.66   Mean   :1.095  
##  3rd Qu.:51.00                      3rd Qu.:34.69   3rd Qu.:2.000  
##  Max.   :64.00                      Max.   :53.13   Max.   :5.000  
##     smoker             region             charges     
##  Length:1338        Length:1338        Min.   : 1122  
##  Class :character   Class :character   1st Qu.: 4740  
##  Mode  :character   Mode  :character   Median : 9382  
##                                        Mean   :13270  
##                                        3rd Qu.:16640  
##                                        Max.   :63770
skim(dbs) # Analisis mas detallado
Data summary
Name dbs
Number of rows 1338
Number of columns 7
_______________________
Column type frequency:
character 3
numeric 4
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
sex 0 1 4 6 0 2 0
smoker 0 1 2 3 0 2 0
region 0 1 9 9 0 4 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 0 1 39.21 14.05 18.00 27.00 39.00 51.00 64.00 ▇▅▅▆▆
bmi 0 1 30.66 6.10 15.96 26.30 30.40 34.69 53.13 ▂▇▇▂▁
children 0 1 1.09 1.21 0.00 0.00 1.00 2.00 5.00 ▇▂▂▁▁
charges 0 1 13270.42 12110.01 1121.87 4740.29 9382.03 16639.91 63770.43 ▇▂▁▁▁

Note que inicialmente no tenemos valores faltantes en ninguna columna, pero cuando hagamos modificaciones ese numero cambiara

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

# Hacemos boxplots para analizar posibilidad de outliers 
dbs %>%
  select(where(is.numeric)) %>%  # Ejemplo excluyendo Outcome
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = "", y = valor)) +
  geom_boxplot(fill = "lightgreen") +
  facet_wrap(~variable, scales = "free_y") +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())

# Barras para variables categóricas
dbs %>%
  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") %>%
  ggplot(aes(x = fct_reorder(valor, n), y = n, fill = variable)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  facet_wrap(~variable, scales = "free") +
  theme_minimal()

# Torta para variables categóricas
dbs %>%
  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") %>%
  ggplot(aes(x = "", y = n, fill = valor)) +
  geom_col(width = 1, color = "white") +
  coord_polar(theta = "y") +
  facet_wrap(~variable) +
  theme_void()

Transformamos los valores pedidos a faltantes

dbs$bmi[dbs$bmi > 40] <- NA
dbs$charges[dbs$charges == 0] <- NA

Revisamos por normalidad

dbs %>%
  select(where(is.numeric)) %>%
  summarise(across(everything(), 
                   ~ shapiro.test(.x)$p.value))
##            age          bmi     children      charges
## 1 5.692046e-22 9.549912e-09 5.066436e-36 1.150523e-36
# Note que se rechaza la normalidad en todas

Ahora tratamos revisamos por datos atipicos

# Empezamos tratando nuestros atipicos usando el test de grubbs en age y children

grubbs.test(dbs$age, type = 10, opposite = FALSE, two.sided = FALSE)
## 
##  Grubbs test for one outlier
## 
## data:  dbs$age
## G = 1.76463, U = 0.99767, p-value = 1
## alternative hypothesis: highest value 64 is an outlier
# No se rechaza la hipotesis nula de que el valor mas alto no es atipico. Volvemos a hacer el test con el mas bajo
grubbs.test(dbs$age, type = 10, opposite = TRUE)
## 
##  Grubbs test for one outlier
## 
## data:  dbs$age
## G = 1.50940, U = 0.99829, p-value = 1
## alternative hypothesis: lowest value 18 is an outlier
# Tampoco se rechaza la hipotesis nula de que el valor mas bajo no es atipico.

grubbs.test(dbs$children, type = 10, opposite = FALSE, two.sided = FALSE)
## 
##  Grubbs test for one outlier
## 
## data:  dbs$children
## G = 3.23941, U = 0.99215, p-value = 0.7852
## alternative hypothesis: highest value 5 is an outlier
# No se rechaza la hipotesis nula de que el valor mas alto no es atipico.
# No sera necesario repetirlo con el mas bajo (0)

# Para age y children no hay atipicos

boxplot.stats(dbs$bmi)$out
## numeric(0)
# Acorde al boxplot bmi tampoco tiene valores atipicos

# Usaremos el filtro de Hamplin para revisar por outliers en charges
lower_bound <- median(dbs$charges, na.rm = TRUE) - 3 * mad(dbs$Insulin, constant = 1, na.rm = TRUE)
upper_bound <- median(dbs$charges, na.rm = TRUE) + 3 * mad(dbs$Insulin, constant = 1, na.rm = TRUE)
which(dbs$charges < lower_bound | dbs$charges > upper_bound)
## integer(0)
# Acorde al filtro de Hamplin, tampoco hay outliers en charges

En conclusion ninguna variable tiene outliers por lo cual no sera necesario tratarlos

gg_miss_var(dbs, show_pct = TRUE)

missmap(dbs, col = c("red", "blue"), legend = TRUE)

Como no hay una sola instancia con charges == 0, el porcentaje de faltantes para este se mantiene en 0 Mientras tanto el de bmi sube a 8% por lo que tendremos que tratalos. Aunque 8% es un porcentaje lo suficientemente bajo como para poder imputar por la mediana, optaremos por usar diversos metodos de mice y comparar cual es el mejor

metodos <- c("pmm", "norm.predict", "norm.nob", "norm")
imputaciones <- list()

# Guardamos la distribución original sin NA
bmi_original <- na.omit(dbs$bmi)

# Histograma original
ggplot(data.frame(bmi = bmi_original), aes(x = bmi)) +
  geom_histogram(binwidth = 1, fill = "red", alpha = 0.6) +
  labs(title = "Distribución original de BMI (sin NA)") +
  theme_minimal()

# Funcion para graficar comparacion directa del metodo de imputacion y la columna sin NA's
graficar_comparacion <- function(original, imputado, variable, metodo) {
  df_plot <- data.frame(
    valor = c(original, imputado),
    tipo = rep(c("Original", paste("Imputado -", metodo)), each = length(original))
  )
  ggplot(df_plot, aes(x = valor, fill = tipo)) +
    geom_density(alpha = 0.4) +
    labs(title = paste("Distribución:", variable),
         subtitle = paste("Método:", metodo)) +
    theme_minimal()
}

plots_bmi <- list()

# Ciclo para aplicar los metodos y guardar graficas
for (m in metodos) {
  # Imputar SOLO bmi
  imp <- mice(dbs, method = m, m = 1, maxit = 5, seed = 123, printFlag = FALSE)
  dbs_imp <- complete(imp)
  imputaciones[[m]] <- dbs_imp
  
  plots_bmi[[m]] <- graficar_comparacion(dbs$bmi, imputaciones[[m]]$bmi, "bmi", m)
  
}
## Warning: Number of logged events: 3
## Warning: Number of logged events: 3
## Warning: Number of logged events: 3
## Warning: Number of logged events: 3
# Observamos las graficas
wrap_plots(plots_bmi)
## Warning: Removed 91 rows containing non-finite outside the scale range
## (`stat_density()`).
## Warning: Removed 91 rows containing non-finite outside the scale range
## (`stat_density()`).
## Removed 91 rows containing non-finite outside the scale range
## (`stat_density()`).
## Removed 91 rows containing non-finite outside the scale range
## (`stat_density()`).

Ahora hacemos pruebas de normalidad de nuevo para definir que pruebas podemos hacer

print(shapiro.test(bmi_original))
## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Shapiro-Wilk normality test
## 
## Data:                            bmi_original
## 
## Test Statistic:                  W = 0.9877055
## 
## P-value:                         9.549912e-09
for (m in metodos) {
  cat("\nMétodo:", m, "\n")
  print(shapiro.test(imputaciones[[m]]$bmi))
}
## 
## Método: pmm 
## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Shapiro-Wilk normality test
## 
## Data:                            imputaciones[[m]]$bmi
## 
## Test Statistic:                  W = 0.9869868
## 
## P-value:                         1.436106e-09
## 
## 
## Método: norm.predict 
## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Shapiro-Wilk normality test
## 
## Data:                            imputaciones[[m]]$bmi
## 
## Test Statistic:                  W = 0.9904011
## 
## P-value:                         1.131518e-07
## 
## 
## Método: norm.nob 
## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Shapiro-Wilk normality test
## 
## Data:                            imputaciones[[m]]$bmi
## 
## Test Statistic:                  W = 0.9889521
## 
## P-value:                         1.614491e-08
## 
## 
## Método: norm 
## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Shapiro-Wilk normality test
## 
## Data:                            imputaciones[[m]]$bmi
## 
## Test Statistic:                  W = 0.9907538
## 
## P-value:                         1.862311e-07

Vemos que todas siguen siendo no normales entonces aplicamos una prueba no parametrica Nuestras hipotesis son:

Hipótesis nula (H0): Las dos distribuciones tienen la misma mediana (o provienen de poblaciones con la misma distribución).

Hipótesis alternativa (H1): Las distribuciones difieren en su mediana (o en su ubicación).

for (m in metodos){
  # Prueba no paramétrica (Wilcoxon)
  test <- wilcox.test(bmi_original, imputaciones[[m]]$bmi)
  cat("\nMétodo:", m, "- p-valor Wilcoxon:", test$p.value, "\n")
}
## 
## Método: pmm - p-valor Wilcoxon: 0.8121035 
## 
## Método: norm.predict - p-valor Wilcoxon: 0.9333836 
## 
## Método: norm.nob - p-valor Wilcoxon: 0.8557098 
## 
## Método: norm - p-valor Wilcoxon: 0.757973

Para ningun metodo se rechaza la hipotesis nula entonces no encontramos evidencia suficiente para decir que la imputación cambió la distribución de bmi en términos de su mediana.

Como todos los metodos son validos, preferiremos pmm para imputar, ya que es la que se ve mas robusta en las graficas