TUTORIAL: IMPUTAÇÃO DE DADOS FALTANTES COM RANDOM FOREST

PASSO 1: CARREGAR PACOTES NECESSÁRIOS

pacotes <- c("tidyverse", "missForest", "naniar", "VIM", "patchwork", "randomForest")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
  }
}
lapply(pacotes, library, character.only = T)

[[1]]
 [1] "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
 [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics" 
[13] "grDevices" "utils"     "datasets"  "methods"   "base"     

[[2]]
 [1] "missForest" "lubridate"  "forcats"    "stringr"    "dplyr"     
 [6] "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"   
[11] "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"     
[16] "datasets"   "methods"    "base"      

[[3]]
 [1] "naniar"     "missForest" "lubridate"  "forcats"    "stringr"   
 [6] "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"    
[11] "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices" 
[16] "utils"      "datasets"   "methods"    "base"      

[[4]]
 [1] "VIM"        "grid"       "colorspace" "naniar"     "missForest"
 [6] "lubridate"  "forcats"    "stringr"    "dplyr"      "purrr"     
[11] "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse" 
[16] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[21] "methods"    "base"      

[[5]]
 [1] "patchwork"  "VIM"        "grid"       "colorspace" "naniar"    
 [6] "missForest" "lubridate"  "forcats"    "stringr"    "dplyr"     
[11] "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"   
[16] "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"     
[21] "datasets"   "methods"    "base"      

[[6]]
 [1] "randomForest" "patchwork"    "VIM"          "grid"         "colorspace"  
 [6] "naniar"       "missForest"   "lubridate"    "forcats"      "stringr"     
[11] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
[16] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
[21] "utils"        "datasets"     "methods"      "base"        

PASSO 2: CRIAR DADOS DE EXEMPLO COM MISSINGS

Vamos criar um dataset que simula dados reais com relações não-lineares entre variáveis.

n <- 500  # Número de observações

dados_completos <- tibble(
  # Variável idade (20 a 70 anos)
  idade = runif(n, 20, 70),
  # Renda depende da idade de forma não-linear (curva quadrática)
  renda = 2000 + 500 * idade - 3 * idade^2 + rnorm(n, 0, 5000),
  # Anos de estudo (relacionado com idade e renda)
  anos_estudo = 8 + 0.15 * idade + 0.0001 * renda + rnorm(n, 0, 2),
  # Pontuação de crédito (depende de várias variáveis)
  score_credito = 300 + 5 * anos_estudo + 0.01 * renda + 
    2 * idade + rnorm(n, 0, 50),
  # Gasto mensal (interação complexa)
  gasto_mensal = 500 + 0.3 * renda - 10 * anos_estudo + 
    0.05 * score_credito + rnorm(n, 0, 1000),
  # Variável categórica: região
  regiao = sample(c("Norte", "Sul", "Leste", "Oeste"), n, replace = TRUE)
) |> 
  # Garantir valores positivos e limites realistas
  mutate(
    renda = pmax(renda, 1000),
    anos_estudo = pmax(pmin(anos_estudo, 20), 0),
    score_credito = pmax(pmin(score_credito, 850), 300),
    gasto_mensal = pmax(gasto_mensal, 0)
  )

# Visualizar os dados completos
glimpse(dados_completos)

Rows: 500
Columns: 6
$ idade         <dbl> 65.70600, 53.25204, 51.95191, 28.42161, 51.79827, 28.589…
$ renda         <dbl> 21210.249, 14151.038, 22750.717, 21889.524, 24653.446, 1…
$ anos_estudo   <dbl> 18.45166, 18.30477, 14.26286, 15.40780, 19.79574, 7.6773…
$ score_credito <dbl> 764.8162, 662.5703, 742.9540, 634.2709, 757.0938, 590.25…
$ gasto_mensal  <dbl> 6381.3801, 4341.0734, 7406.1504, 5843.9599, 6717.9034, 4…
$ regiao        <chr> "Sul", "Sul", "Leste", "Norte", "Sul", "Sul", "Oeste", "…

PASSO 3: INTRODUZIR DADOS FALTANTES DE FORMA REALISTA

Criar cópia dos dados e introduzir missings usando MCAR (Missing Completely At Random) em algumas variáveis e MAR (Missing At Random) em outras - mais realista.

dados_com_missing <- dados_completos |> 
  mutate(
    # 20% de missings em renda (MCAR)
    renda = if_else(runif(n) < 0.20, NA_real_, renda),
    # Missings em anos_estudo dependem da idade (MAR)
    # Pessoas mais jovens tendem a não responder
    anos_estudo = if_else(idade < 30 & runif(n) < 0.25, 
                          NA_real_, anos_estudo),
    # 15% de missings em score_credito (MCAR)
    score_credito = if_else(runif(n) < 0.15, NA_real_, score_credito),
    # Missings em gasto_mensal relacionados com renda alta (MAR)
    gasto_mensal = if_else(renda > 40000 & runif(n) < 0.30 | runif(n) < 0.10,
                           NA_real_, gasto_mensal),
    # 10% de missings em região
    regiao = if_else(runif(n) < 0.10, NA_character_, regiao)
  )

PASSO 4: ANÁLISE EXPLORATÓRIA DOS DADOS FALTANTES

# Contar missings por variável
resumo_missing <- dados_com_missing |> 
  summarise(across(everything(), ~sum(is.na(.)))) |> 
  pivot_longer(everything(), 
               names_to = "variavel", 
               values_to = "n_missing") |> 
  mutate(
    perc_missing = round(n_missing / n * 100, 2),
    total_obs = n
  ) |> 
  arrange(desc(n_missing))

print(resumo_missing)

# A tibble: 6 × 4
  variavel      n_missing perc_missing total_obs
  <chr>             <int>        <dbl>     <dbl>
1 renda               102         20.4       500
2 score_credito        83         16.6       500
3 gasto_mensal         65         13         500
4 regiao               50         10         500
5 anos_estudo          34          6.8       500
6 idade                 0          0         500

# Visualização 1: Padrão de dados faltantes
vis_missing_1 <- ggplot(dados_com_missing, 
                        aes(x = 1:nrow(dados_com_missing))) +
  geom_miss_point(aes(y = renda), alpha = 0.4, size = 2) +
  labs(title = "Padrão de Dados Faltantes: Renda",
       subtitle = "Pontos vermelhos indicam valores ausentes",
       x = "Índice da Observação",
       y = "Renda") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

# Visualização 2: Proporção de missings por variável
vis_missing_2 <- resumo_missing |> 
  filter(n_missing > 0) |> 
  ggplot(aes(x = reorder(variavel, perc_missing), 
             y = perc_missing,
             fill = perc_missing)) +
  geom_col() +
  geom_text(aes(label = paste0(perc_missing, "%")), 
            hjust = -0.2, size = 3.5) +
  coord_flip() +
  scale_fill_gradient(low = "#FEE5D9", high = "#A50F15") +
  labs(title = "Proporção de Dados Faltantes por Variável",
       x = NULL,
       y = "Percentual de Valores Ausentes (%)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold")) +
  ylim(0, max(resumo_missing$perc_missing) * 1.15)

# Visualização 3: Matriz de padrões de missings
vis_missing_3 <- gg_miss_upset(dados_com_missing, 
                               nsets = 6,
                               nintersects = 15)

# Combinar visualizações
print(vis_missing_1)

print(vis_missing_2)

print(vis_missing_3)

PASSO 5: IMPUTAÇÃO COM RANDOM FOREST

# Converter região para factor (necessário para missForest)
dados_para_imputar <- dados_com_missing |> 
  mutate(regiao = as.factor(regiao))

Para executar imputação com missForest estabeleço os seguintes parâmetros: - ntree: número de árvores na floresta
- mtry: número de variáveis testadas em cada divisão
- maxiter: máximo de iterações
- verbose: mostrar progresso

resultado_imputacao <- missForest(
  xmis = as.data.frame(dados_para_imputar),
  ntree = 100,           # 100 árvores para cada floresta
  mtry = 2,              # Raiz quadrada do número de variáveis
  maxiter = 10,          # Máximo de 10 iterações
  verbose = TRUE         # Mostrar progresso
)

  missForest iteration 1 in progress...done!
    estimated error(s): 0.1802504 0.72 
    difference(s): 0.01215697 0.074 
    time: 0.16 seconds

  missForest iteration 2 in progress...done!
    estimated error(s): 0.176472 0.6866667 
    difference(s): 0.0002999954 0.046 
    time: 0.14 seconds

  missForest iteration 3 in progress...done!
    estimated error(s): 0.1747374 0.7088889 
    difference(s): 0.0002331696 0.036 
    time: 0.14 seconds

  missForest iteration 4 in progress...done!
    estimated error(s): 0.1769141 0.6933333 
    difference(s): 0.0001003542 0.026 
    time: 0.12 seconds

  missForest iteration 5 in progress...done!
    estimated error(s): 0.1778907 0.7 
    difference(s): 0.0001084176 0.018 
    time: 0.14 seconds

  missForest iteration 6 in progress...done!
    estimated error(s): 0.1781782 0.7022222 
    difference(s): 8.033496e-05 0.018 
    time: 0.14 seconds

  missForest iteration 7 in progress...done!
    estimated error(s): 0.1745193 0.7088889 
    difference(s): 9.338581e-05 0.018 
    time: 0.13 seconds

Agora vamos extrair dados imputados

dados_imputados <- as_tibble(resultado_imputacao$ximp) |> 
  mutate(regiao = as.character(regiao))

Informações de Erros de imputação

cat("NRMSE (Normalized Root Mean Squared Error) para variáveis numéricas:\n")

NRMSE (Normalized Root Mean Squared Error) para variáveis numéricas:

print(resultado_imputacao$OOBerror[1])  # Quanto menor, melhor

    NRMSE 
0.1781782 

cat("\nPFC (Proportion of Falsely Classified) para variáveis categóricas:\n")

PFC (Proportion of Falsely Classified) para variáveis categóricas:

print(resultado_imputacao$OOBerror[2])  # Quanto menor, melhor

      PFC 
0.7022222 

PASSO 6: VALIDAÇÃO DA IMPUTAÇÃO

Vamos comparar a distribuição das variáveis antes e depois da imputação.

# Identificar quais valores foram imputados
dados_comparacao <- dados_completos |> 
  mutate(origem = "Original") |> 
  bind_rows(
    dados_imputados |> 
      mutate(origem = "Imputado")
  ) |> 
  mutate(id = rep(1:n, 2)) |> 
  # Criar flag para valores que foram imputados
  left_join(
    dados_com_missing |> 
      mutate(id = row_number()) |> 
      reframe(
        id = id,
        renda_era_missing = is.na(renda),
        anos_estudo_era_missing = is.na(anos_estudo),
        score_era_missing = is.na(score_credito),
        gasto_era_missing = is.na(gasto_mensal)
      ),
    by = "id"
  )


# Visualização 4: Comparar distribuições - Renda
vis_validacao_1 <- dados_comparacao |> 
  filter(origem == "Original" | (origem == "Imputado" & renda_era_missing)) |> 
  ggplot(aes(x = renda, fill = origem)) +
  geom_density(alpha = 0.6) +
  scale_fill_manual(values = c("Original" = "#2166AC", 
                               "Imputado" = "#B2182B")) +
  labs(title = "Distribuição da Renda: Original vs Imputado",
       subtitle = "Os valores imputados seguem a distribuição original?",
       x = "Renda",
       y = "Densidade",
       fill = "Tipo") +
  theme_classic(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = "bottom")

# Visualização 5: Comparar distribuições - Score de Crédito
vis_validacao_2 <- dados_comparacao |> 
  filter(origem == "Original" | (origem == "Imputado" & score_era_missing)) |> 
  ggplot(aes(x = score_credito, fill = origem)) +
  geom_density(alpha = 0.6) +
  scale_fill_manual(values = c("Original" = "#2166AC", 
                               "Imputado" = "#B2182B")) +
  labs(title = "Distribuição do Score: Original vs Imputado",
       x = "Score de Crédito",
       y = "Densidade",
       fill = "Tipo") +
  theme_classic(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = "bottom")

# Visualização 6: Boxplots comparativos
vis_validacao_3 <- dados_comparacao |> 
  filter(origem == "Original" | (origem == "Imputado" & gasto_era_missing)) |> 
  ggplot(aes(x = origem, y = gasto_mensal, fill = origem)) +
  geom_boxplot(alpha = 0.7) +
  scale_fill_manual(values = c("Original" = "#2166AC", 
                               "Imputado" = "#B2182B")) +
  labs(title = "Gasto Mensal: Original vs Imputado",
       x = NULL,
       y = "Gasto Mensal") +
  theme_classic(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = "none")

# Visualização 7: Scatter plot - Relação entre variáveis
vis_validacao_4 <- dados_comparacao |> 
  filter(origem == "Imputado") |> 
  ggplot(aes(x = idade, y = renda)) +
  geom_point(aes(color = renda_era_missing), alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "black", 
              linewidth = 1, linetype = "dashed") +
  scale_color_manual(
    values = c("TRUE" = "#B2182B", "FALSE" = "#2166AC"),
    labels = c("TRUE" = "Valor Imputado", "FALSE" = "Valor Original")
  ) +
  labs(title = "Relação Idade vs Renda",
       subtitle = "Valores imputados preservam a relação não-linear",
       x = "Idade",
       y = "Renda",
       color = NULL) +
  theme_classic(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = "bottom")

# Imprimir visualizações de validação
print(vis_validacao_1)

print(vis_validacao_2)

print(vis_validacao_3)

print(vis_validacao_4)

`geom_smooth()` using formula = 'y ~ x'

PASSO 7: MÉTRICAS QUANTITATIVAS DE VALIDAÇÃO

# Função para calcular RMSE
calcular_rmse <- function(original, imputado, indices_missing) {
  valores_originais <- original[indices_missing]
  valores_imputados <- imputado[indices_missing]
  sqrt(mean((valores_originais - valores_imputados)^2))
}

# Função para calcular MAE
calcular_mae <- function(original, imputado, indices_missing) {
  valores_originais <- original[indices_missing]
  valores_imputados <- imputado[indices_missing]
  mean(abs(valores_originais - valores_imputados))
}

# Calcular métricas para cada variável
indices_missing_renda <- which(is.na(dados_com_missing$renda))
indices_missing_anos <- which(is.na(dados_com_missing$anos_estudo))
indices_missing_score <- which(is.na(dados_com_missing$score_credito))
indices_missing_gasto <- which(is.na(dados_com_missing$gasto_mensal))

metricas_validacao <- tibble(
  variavel = c("Renda", "Anos de Estudo", "Score de Crédito", "Gasto Mensal"),
  rmse = c(
    calcular_rmse(dados_completos$renda, dados_imputados$renda, 
                  indices_missing_renda),
    calcular_rmse(dados_completos$anos_estudo, dados_imputados$anos_estudo, 
                  indices_missing_anos),
    calcular_rmse(dados_completos$score_credito, dados_imputados$score_credito, 
                  indices_missing_score),
    calcular_rmse(dados_completos$gasto_mensal, dados_imputados$gasto_mensal, 
                  indices_missing_gasto)
  ),
  mae = c(
    calcular_mae(dados_completos$renda, dados_imputados$renda, 
                 indices_missing_renda),
    calcular_mae(dados_completos$anos_estudo, dados_imputados$anos_estudo, 
                 indices_missing_anos),
    calcular_mae(dados_completos$score_credito, dados_imputados$score_credito, 
                 indices_missing_score),
    calcular_mae(dados_completos$gasto_mensal, dados_imputados$gasto_mensal, 
                 indices_missing_gasto)
  )
) %>%
  mutate(
    rmse = round(rmse, 2),
    mae = round(mae, 2)
  )

print(metricas_validacao)

# A tibble: 4 × 3
  variavel            rmse     mae
  <chr>              <dbl>   <dbl>
1 Renda            3126.   2445.  
2 Anos de Estudo      2.27    1.79
3 Score de Crédito   56.5    47.3 
4 Gasto Mensal     1229.    920.  

# Visualização 8: Métricas de erro
vis_metricas <- metricas_validacao %>%
  pivot_longer(cols = c(rmse, mae), 
               names_to = "metrica", 
               values_to = "valor") %>%
  ggplot(aes(x = variavel, y = valor, fill = metrica)) +
  geom_col(position = "dodge", alpha = 0.8) +
  geom_text(aes(label = round(valor, 1)), 
            position = position_dodge(width = 0.9),
            vjust = -0.5, size = 3) +
  scale_fill_manual(
    values = c("rmse" = "#D6604D", "mae" = "#4393C3"),
    labels = c("rmse" = "RMSE", "mae" = "MAE")
  ) +
  labs(title = "Erro de Imputação por Variável",
       subtitle = "RMSE = Root Mean Squared Error | MAE = Mean Absolute Error",
       x = NULL,
       y = "Valor do Erro",
       fill = "Métrica") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

print(vis_metricas)

PASSO 8: COMPARAÇÃO COM MÉTODO SIMPLES (MÉDIA)

# Imputação simples por média
dados_imputados_media <- dados_com_missing %>%
  mutate(
    renda = if_else(is.na(renda), mean(renda, na.rm = TRUE), renda),
    anos_estudo = if_else(is.na(anos_estudo), 
                          mean(anos_estudo, na.rm = TRUE), anos_estudo),
    score_credito = if_else(is.na(score_credito), 
                            mean(score_credito, na.rm = TRUE), score_credito),
    gasto_mensal = if_else(is.na(gasto_mensal), 
                           mean(gasto_mensal, na.rm = TRUE), gasto_mensal)
  )

# Calcular métricas para imputação por média
metricas_media <- tibble(
  variavel = c("Renda", "Anos de Estudo", "Score de Crédito", "Gasto Mensal"),
  rmse_rf = metricas_validacao$rmse,
  rmse_media = c(
    calcular_rmse(dados_completos$renda, dados_imputados_media$renda, 
                  indices_missing_renda),
    calcular_rmse(dados_completos$anos_estudo, 
                  dados_imputados_media$anos_estudo, indices_missing_anos),
    calcular_rmse(dados_completos$score_credito, 
                  dados_imputados_media$score_credito, indices_missing_score),
    calcular_rmse(dados_completos$gasto_mensal, 
                  dados_imputados_media$gasto_mensal, indices_missing_gasto)
  )
) %>%
  mutate(
    rmse_media = round(rmse_media, 2),
    melhoria_percentual = round((rmse_media - rmse_rf) / rmse_media * 100, 1)
  )

print(metricas_media)

# A tibble: 4 × 4
  variavel         rmse_rf rmse_media melhoria_percentual
  <chr>              <dbl>      <dbl>               <dbl>
1 Renda            3126.      6335.                  50.6
2 Anos de Estudo      2.27       4.38                48.2
3 Score de Crédito   56.5      100.                  43.5
4 Gasto Mensal     1229.      2047.                  40  

# Visualização 9: Comparação de métodos
vis_comparacao <- metricas_media %>%
  select(variavel, rmse_rf, rmse_media) %>%
  pivot_longer(cols = starts_with("rmse"), 
               names_to = "metodo", 
               values_to = "rmse") %>%
  mutate(metodo = recode(metodo,
                         "rmse_rf" = "Random Forest",
                         "rmse_media" = "Média Simples")) %>%
  ggplot(aes(x = variavel, y = rmse, fill = metodo)) +
  geom_col(position = "dodge", alpha = 0.8) +
  geom_text(aes(label = round(rmse, 1)), 
            position = position_dodge(width = 0.9),
            vjust = -0.5, size = 3) +
  scale_fill_manual(values = c("Random Forest" = "#1B9E77", 
                               "Média Simples" = "#D95F02")) +
  labs(title = "Comparação de Métodos de Imputação",
       subtitle = "Random Forest consistentemente supera imputação por média",
       x = NULL,
       y = "RMSE",
       fill = "Método") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

print(vis_comparacao)

RESUMO DO PROCESSO DE IMPUTAÇÃO

cat("✓ Dados originais criados:", n, "observações\n")

✓ Dados originais criados: 500 observações

cat("✓ Variáveis com missings introduzidos:", 
    sum(resumo_missing$n_missing > 0), "\n")

✓ Variáveis com missings introduzidos: 5 

cat("✓ Total de valores faltantes:", sum(resumo_missing$n_missing), "\n")

✓ Total de valores faltantes: 334 

cat("✓ Método de imputação: Random Forest (missForest)\n")

✓ Método de imputação: Random Forest (missForest)

cat("✓ Número de árvores utilizadas: 100\n")

✓ Número de árvores utilizadas: 100

cat("✓ NRMSE médio:", round(resultado_imputacao$OOBerror[1], 4), "\n")

✓ NRMSE médio: 0.1782 

cat("✓ Melhoria média sobre imputação por média:", 
    round(mean(metricas_media$melhoria_percentual), 1), "%\n\n")

✓ Melhoria média sobre imputação por média: 45.6 %

CONCLUSÃO

Random Forest capturou com sucesso as relações não-lineares presentes nos dados, gerando imputações que preservam a estrutura e distribuição dos valores originais. Este método é superior à imputação simples por média, especialmente quando há interações complexas entre variáveis.

PASSO 9: EXPORTAR DADOS IMPUTADOS

# Salvar dados imputados
write_csv(dados_imputados, "dados_imputados_rf.csv")