Previsão de AVC
Arnaldo Chebl
2023-11-20
Descrição do problema e análise das métricas de desempenho
O objetivo do trabalho é prever se o paciente vai ter AVC, antes que ele de fato o tenha. Vamos utilizar o framework ‘tidymodels’ para modelar o dataset, empregando técnicas de ‘Boosting’ e ‘Random Forest’, ambas com e sem ajuste de parâmetros (‘Tuning’). Além disso, usaremos o Keras para desenvolver um modelo baseado em rede neural. Nosso foco principal será minimizar falsos negativos, visando identificar com antecedência os pacientes em risco de derrame para prevenir fatalidades.
Setup e Bibliotecas
Clique para expandir e ver as bibliotecas carregadas
library(tidymodels)
library(skimr)
library(factoextra)
library(ggplot2)
library(plotly)
library(dplyr)
library(doParallel)
library(vip)
library(broom)
library(knitr)
library("reticulate")
reticulate::use_python("C:/Users/arnch/anaconda3/python.exe")
library(keras)
mnist <- dataset_mnist()
library(pROC)
library(conflicted)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")Input da Base e ajustes básicos
Verifica tipos e nulos
Clique para expandir
df <- read.csv("healthcare-dataset-stroke-data.csv")
df %>%
lapply(type_sum) %>%
as_tibble() %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Tipo") %>%
inner_join(
df %>%
summarise(across(everything(), ~sum(is.na(.)))) %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Total NA")
)## # A tibble: 12 × 3
## Coluna Tipo `Total NA`
## <chr> <chr> <int>
## 1 id int 0
## 2 gender chr 0
## 3 age dbl 0
## 4 hypertension int 0
## 5 heart_disease int 0
## 6 ever_married chr 0
## 7 work_type chr 0
## 8 Residence_type chr 0
## 9 avg_glucose_level dbl 0
## 10 bmi chr 0
## 11 smoking_status chr 0
## 12 stroke int 0df %>% glimpse()## Rows: 5,110
## Columns: 12
## $ id <int> 9046, 51676, 31112, 60182, 1665, 56669, 53882, 10434…
## $ gender <chr> "Male", "Female", "Male", "Female", "Female", "Male"…
## $ age <dbl> 67, 61, 80, 49, 79, 81, 74, 69, 59, 78, 81, 61, 54, …
## $ hypertension <int> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1…
## $ heart_disease <int> 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0…
## $ ever_married <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "No…
## $ work_type <chr> "Private", "Self-employed", "Private", "Private", "S…
## $ Residence_type <chr> "Urban", "Rural", "Rural", "Urban", "Rural", "Urban"…
## $ avg_glucose_level <dbl> 228.69, 202.21, 105.92, 171.23, 174.12, 186.21, 70.0…
## $ bmi <chr> "36.6", "N/A", "32.5", "34.4", "24", "29", "27.4", "…
## $ smoking_status <chr> "formerly smoked", "never smoked", "never smoked", "…
## $ stroke <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…df %>% skim() %>% View()
print(paste("Percentual de casos com derrames:", 100 * (df %>% filter(stroke == 1) %>% nrow() / nrow(df))))## [1] "Percentual de casos com derrames: 4.87279843444227"Arruma os tipos das variáveis
Clique para expandir
df <- df %>%
mutate(
gender = as.factor(gender),
hypertension = as.factor(hypertension),
heart_disease = as.factor(heart_disease),
ever_married = as.factor(ever_married),
work_type = as.factor(work_type),
Residence_type = as.factor(Residence_type),
smoking_status = as.factor(smoking_status),
stroke = as.factor(stroke)
)
df |> glimpse()## Rows: 5,110
## Columns: 12
## $ id <int> 9046, 51676, 31112, 60182, 1665, 56669, 53882, 10434…
## $ gender <fct> Male, Female, Male, Female, Female, Male, Male, Fema…
## $ age <dbl> 67, 61, 80, 49, 79, 81, 74, 69, 59, 78, 81, 61, 54, …
## $ hypertension <fct> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1…
## $ heart_disease <fct> 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0…
## $ ever_married <fct> Yes, Yes, Yes, Yes, Yes, Yes, Yes, No, Yes, Yes, Yes…
## $ work_type <fct> Private, Self-employed, Private, Private, Self-emplo…
## $ Residence_type <fct> Urban, Rural, Rural, Urban, Rural, Urban, Rural, Urb…
## $ avg_glucose_level <dbl> 228.69, 202.21, 105.92, 171.23, 174.12, 186.21, 70.0…
## $ bmi <chr> "36.6", "N/A", "32.5", "34.4", "24", "29", "27.4", "…
## $ smoking_status <fct> formerly smoked, never smoked, never smoked, smokes,…
## $ stroke <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…#Abaixo verificamos que a variável 'bmi', que está como character, possui valores 'N/A'.
#Vamos criar uma flag para captar a informação dos valores missing desta váriavel e vamos transformar BMI no tipo numérico.
df$bmi %>% table() %>% sort(decreasing = TRUE) %>% head()## .
## N/A 28.7 28.4 26.1 26.7 27.6
## 201 41 38 37 37 37#criamos feature flag
df$flag <- ifelse(df$bmi == "N/A", 1, 0)
df$bmi <- ifelse(df$bmi == "N/A", NA, df$bmi)
df <- df %>%
mutate(
flag = as.factor(flag),
bmi = as.double(bmi),
)
df %>% head(5) %>% glimpse()## Rows: 5
## Columns: 13
## $ id <int> 9046, 51676, 31112, 60182, 1665
## $ gender <fct> Male, Female, Male, Female, Female
## $ age <dbl> 67, 61, 80, 49, 79
## $ hypertension <fct> 0, 0, 0, 0, 1
## $ heart_disease <fct> 1, 0, 1, 0, 0
## $ ever_married <fct> Yes, Yes, Yes, Yes, Yes
## $ work_type <fct> Private, Self-employed, Private, Private, Self-emplo…
## $ Residence_type <fct> Urban, Rural, Rural, Urban, Rural
## $ avg_glucose_level <dbl> 228.7, 202.2, 105.9, 171.2, 174.1
## $ bmi <dbl> 36.6, NA, 32.5, 34.4, 24.0
## $ smoking_status <fct> formerly smoked, never smoked, never smoked, smokes,…
## $ stroke <fct> 1, 1, 1, 1, 1
## $ flag <fct> 0, 1, 0, 0, 0Análise das categorias
# input de dados ----------------------------------------------------------
categorical_variables <- df %>%
select_if(is.factor) %>%
names()
for (i in categorical_variables) {
cat("Variable:", i, "\n")
print(table(df[[i]]))
cat("\n")
}## Variable: gender
##
## Female Male Other
## 2994 2115 1
##
## Variable: hypertension
##
## 0 1
## 4612 498
##
## Variable: heart_disease
##
## 0 1
## 4834 276
##
## Variable: ever_married
##
## No Yes
## 1757 3353
##
## Variable: work_type
##
## children Govt_job Never_worked Private Self-employed
## 687 657 22 2925 819
##
## Variable: Residence_type
##
## Rural Urban
## 2514 2596
##
## Variable: smoking_status
##
## formerly smoked never smoked smokes Unknown
## 885 1892 789 1544
##
## Variable: stroke
##
## 0 1
## 4861 249
##
## Variable: flag
##
## 0 1
## 4909 201#há apenas um genero do tipo "outro", iremos excluí-lo pois não será possivel treinar o modelo com apenas um registro
df <- df %>% plotly::filter(gender != "Other")Tratando nulos da variável ‘bmi’ com a mediana
df %>%
lapply(type_sum) %>%
as_tibble() %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Tipo") %>%
inner_join(
df %>%
summarise(across(everything(), ~sum(is.na(.)))) %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Total NA")
)## # A tibble: 13 × 3
## Coluna Tipo `Total NA`
## <chr> <chr> <int>
## 1 id int 0
## 2 gender fct 0
## 3 age dbl 0
## 4 hypertension fct 0
## 5 heart_disease fct 0
## 6 ever_married fct 0
## 7 work_type fct 0
## 8 Residence_type fct 0
## 9 avg_glucose_level dbl 0
## 10 bmi dbl 201
## 11 smoking_status fct 0
## 12 stroke fct 0
## 13 flag fct 0# Aplicar mediana aos nulos
df$bmi[is.na(df$bmi)] <- median(df$bmi, na.rm=TRUE)
# Verificar nulos novamente
df %>%
lapply(type_sum) %>%
as_tibble() %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Tipo") %>%
inner_join(
df %>%
summarise(across(everything(), ~sum(is.na(.)))) %>%
pivot_longer(cols = 1:ncol(df),
names_to = "Coluna",
values_to = "Total NA")
)## # A tibble: 13 × 3
## Coluna Tipo `Total NA`
## <chr> <chr> <int>
## 1 id int 0
## 2 gender fct 0
## 3 age dbl 0
## 4 hypertension fct 0
## 5 heart_disease fct 0
## 6 ever_married fct 0
## 7 work_type fct 0
## 8 Residence_type fct 0
## 9 avg_glucose_level dbl 0
## 10 bmi dbl 0
## 11 smoking_status fct 0
## 12 stroke fct 0
## 13 flag fct 0Clusterização: k-means
Análise dos clusters
## # A tibble: 4 × 2
## # Groups: cluster [2]
## cluster pct_stroke
## <fct> <dbl>
## 1 1 0.920
## 2 1 0.0800
## 3 2 0.996
## 4 2 0.00381
Podemos verificar que, em relação ao Cluster 2, no Cluster 1:
- As pessoas são bem mais velhas (a média é 58 x 22)
- As pessoas tem níveis de glicose mais altos (a média é 117 x 91)
- As pessoas tem níveis IMC (BMI) mais altos (a média é 32 x 24)
- AS pessoas sofrem mais infartos (8% x 0.4%)
Análise Exploratória
## [1] "verificamos necessidade de fazer log em 'bmi' e em 'avg_glucose_level', para normaliza-las"
Podemos verificar que:
- Age x BMI: Pessoas no cluster 1 são majoritariamente mais velhas e com valores mais altos de IMC, contendo quase todos os indivíduos que tiveram derrame;
- BMI x avg_glucose_level e age x avg_glucose_level: maior parte dos indivíduos do cluster 2 tem nível médio de glicose abaixo de 150, enquanto no cluster 1 parece haver uma separação mais equilibrada, tendo uma concentração entre 50 e 125, e outra entre 200 e 250;
- Gender: A proporção de homens e mulheres dentro dos grupos de pacientes que tiveram e não tiveram derrame é parecida, indicando que o gênero não possui muita relevância para ocorrência do evento;
- hypertension: Pacientes com hipertensão tem maior proporção dentro do grupo que apresentou derrame em relação ao que não apresentou (26% x 8%);
- heart_disease: Pacientes com doenças cardiacas tem maior proporção dentro do grupo que apresentou derrame em relação ao que não apresentou (18% x 4%);
- ever_married: Pacientes que são ou já foram casados tem maior proporção dentro do grupo que apresentou derrame em relação ao que não apresentou (88% x 64%);
- work_type: Pacientes que trabalham na área privada tem proporções parecidas dentro dos grupos que tiveram e não tiveram derrame, porém autônomos (‘Self-employed’) representam 26% do grupo com derrame enquanto são apenas 15% do grupo que não apresentou;
- Residence_type: A proporção de pacientes que vivem em area urbana e rural dentro dos grupos de que tiveram e não tiveram derrame é parecida, indicando que o o tipo de residência não possui muita relevância para ocorrência do evento;
- smoking_status: formerly smoked representam 28% do grupo com derrame enquanto são apenas 16% do grupo que não apresentou, demais categorias conhecidas possuem proporções semelhantes;
Ajuste dos modelos
Split da base em treino e teste
Clique para expandir
set.seed(234)
split <- initial_split(df, prop = 0.8, strata = stroke)
treinamento <- training(split)
teste <- testing(split)
prop.table(table(treinamento$stroke))##
## 0 1
## 0.95229 0.04771prop.table(table(teste$stroke))##
## 0 1
## 0.94716 0.05284Criação do Tibble de Comparativo do desempenho dos modelos
Clique para expandir
desempenho <- tibble()Criação do recipe, prep e bake do tidymodels com interação
Clique para expandir
# Transformação de dados - formato tidymodels -------------------------------------
set.seed(234)
#faz a receita e prepara
(receita <- recipe(stroke ~ ., data = treinamento) %>%
step_rm(id) %>% #remove colunas selecionadas
step_zv(all_predictors()) %>% #remove variáveis que contém um único valor
step_log(avg_glucose_level, bmi, offset = 0) |> #melhora regressão linear e logistica se variáveis tiverem distribuição exponencial
step_impute_median(bmi) %>% #substitui com mediana casos nulos
step_interact(~ all_predictors():all_predictors()) %>% #cria interações entre todas as variáveis
step_zv(all_predictors()) %>% #remove variáveis que contém um único valor
#step_interact(~ all_predictors():all_predictors()) %>%
step_normalize(all_numeric(), -all_outcomes()) %>% #normaliza
step_other(all_nominal(), -all_outcomes(), threshold = .05, other = "outros") %>% #chama de outros categorias com representatividade abaixo de 5%
step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE) %>% #transforma categoricas em dummies contendo todas categorias nas colunas (one hot)
step_zv(all_predictors()) #remove colunas que contem mesmo valor
)
(receita_prep <- prep(receita))
# obtem os dados de treinamento processados
treinamento_proc <- bake(receita_prep, new_data = NULL)
# obtem os dados de teste processados
teste_proc <- bake(receita_prep, new_data = teste)
base_full_proc <- bake(receita_prep, new_data = df)Regressão Linear - com interação
logreg_cls_spec <-
logistic_reg() %>%
set_engine("glm")
set.seed(1)
fit_lm <- logreg_cls_spec %>% fit(stroke ~ ., data = treinamento_proc)
desempenho <- desempenho %>%
bind_rows(tibble(prob = predict(fit_lm, new_data = teste_proc, type = "prob")$.pred_1,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "LM - interact"))Os p.valores mais significativos de cada combinação de parâmetro e seus betas estimados:
tidy(fit_lm) |> arrange(p.value) |> filter(p.value < 0.05)## # A tibble: 9 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 hypertension1_x_flag1 -0.204 0.0745 -2.74 0.00616
## 2 heart_disease1_x_smoking_statussmokes 0.228 0.0864 2.64 0.00824
## 3 `heart_disease1_x_work_typeSelf-employed` -2.28 0.911 -2.51 0.0122
## 4 heart_disease1_x_work_typePrivate -3.25 1.36 -2.39 0.0169
## 5 heart_disease1_x_work_typeGovt_job -1.57 0.668 -2.35 0.0187
## 6 age_x_smoking_statussmokes -1.19 0.528 -2.26 0.0241
## 7 heart_disease1_x_bmi 2.88 1.35 2.14 0.0324
## 8 smoking_statussmokes_x_flag1 -0.204 0.0991 -2.06 0.0398
## 9 avg_glucose_level_x_flag1 -1.29 0.630 -2.04 0.0410LASSO - com interação
Impacto da regularização na AUC
# Definir os hiperparâmetros para tunagem
LASSO <- logistic_reg(penalty = tune(), mixture = 1) %>%
set_engine("glmnet") %>%
set_mode("classification")
# CV
set.seed(234)
cv_split <- vfold_cv(treinamento, v = 10)
# Paralelização
doParallel::registerDoParallel(makeCluster(16))
tempo <- system.time({
LASSO_grid <- tune_grid(LASSO,
receita,
resamples = cv_split,
grid = 100,
metrics = metric_set(roc_auc))
})
autoplot(LASSO_grid) # Plota os resultados do grid search
Grid dos melhores resultados de acordo com tunning de lambda/penalidade
## # A tibble: 6 × 7
## penalty .metric .estimator mean n std_err .config
## <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.00536 roc_auc binary 0.863 10 0.0118 Preprocessor1_Model078
## 2 0.00425 roc_auc binary 0.863 10 0.0118 Preprocessor1_Model077
## 3 0.00382 roc_auc binary 0.863 10 0.0119 Preprocessor1_Model076
## 4 0.0119 roc_auc binary 0.862 10 0.0115 Preprocessor1_Model081
## 5 0.00266 roc_auc binary 0.862 10 0.0120 Preprocessor1_Model075
## 6 0.00792 roc_auc binary 0.862 10 0.0115 Preprocessor1_Model079Variáveis e interações mais importantes pela biblioteca VIP (que mais afetam o erro do modelo caso sejam alteradas aleatoriamente)
# Seleciona o melhor conjunto de hiperparâmetros
best_LASSO <- LASSO_grid %>%
select_best("roc_auc") #select_by_one_std_err() ou select_best()
#tempo # mostra o tempo de execução
# Faz o fit do modelo com o melhor conjunto de hiperparâmetros
LASSO_fit <- finalize_model(LASSO, parameters = best_LASSO) %>%
fit(stroke ~ ., data = treinamento_proc)
pred_LASSO <- predict(LASSO_fit, new_data = teste_proc, type = "prob")$.pred_1 #pega coluna que dá probabilidade de ter o derrame
desempenho <- desempenho %>%
bind_rows(tibble(prob = pred_LASSO,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "LASSO - interact"))
#desempenho %>% view()
# visualizando variáveis mais importantes
vip(LASSO_fit, num_features = 10, nsim = 10)
#vip::vi(LASSO_fit, num_features = 10, nsim = 10)Interações de variáveis mais importantes: de maior beta em módulo
# Obter os betas em ordem decrescente
tidy(LASSO_fit) |> arrange(desc(abs(estimate))) |> head(8)## # A tibble: 8 × 3
## term estimate penalty
## <chr> <dbl> <dbl>
## 1 (Intercept) -3.79 0.00536
## 2 age_x_avg_glucose_level 1.04 0.00536
## 3 age_x_flag1 0.223 0.00536
## 4 age 0.186 0.00536
## 5 heart_disease1_x_smoking_statussmokes 0.106 0.00536
## 6 smoking_statusUnknown_x_flag1 0.0979 0.00536
## 7 age_x_hypertension1 0.0901 0.00536
## 8 hypertension1_x_flag1 -0.0657 0.00536Criação do recipe, prep e bake do tidymodels sem interação
Clique para expandir
# Transformação de dados - formato tidymodels -------------------------------------
set.seed(234)
#faz a receita e prepara
(receita <- recipe(stroke ~ ., data = treinamento) %>%
step_rm(id) %>% #remove colunas selecionadas
step_zv(all_predictors()) %>% #remove variáveis que contém um único valor
step_log(avg_glucose_level, bmi) |> #melhora regressão linear e logistica se variáveis tiverem crescimento exponencial
step_impute_median(bmi) %>% #substitui com mediana casos nulos
#step_interact(~ all_predictors():all_predictors()) %>% #cria interações entre todas as variáveis
#step_zv(all_predictors()) %>% #remove variáveis que contém um único valor
#step_interact(~ all_predictors():all_predictors()) %>%
step_normalize(all_numeric(), -all_outcomes()) %>% #normaliza
step_other(all_nominal(), -all_outcomes(), threshold = .05, other = "outros") %>% #chama de outros categorias com representatividade abaixo de 5%
step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE) %>% #transforma categoricas em dummies contendo todas categorias nas colunas (one hot)
step_zv(all_predictors()) #remove colunas que contem mesmo valor
)
(receita_prep <- prep(receita))
# obtem os dados de treinamento processados
treinamento_proc <- bake(receita_prep, new_data = NULL)
# obtem os dados de teste processados
teste_proc <- bake(receita_prep, new_data = teste)
base_full_proc <- bake(receita_prep, new_data = df)
#treinamento_proc %>% ncol()
teste %>% ncol()## [1] 14teste_proc %>% ncol()## [1] 27Regressão Linear - sem interação
Os p.valores mais significativos de cada parâmetro e seus betas estimados:
logreg_cls_spec <-
logistic_reg() %>%
set_engine("glm")
set.seed(1)
fit_lm <- logreg_cls_spec %>% fit(stroke ~ ., data = treinamento_proc)
pred_lm <- predict(fit_lm, new_data = teste_proc, type = "prob")$.pred_1
desempenho <- desempenho %>%
bind_rows(tibble(prob = pred_lm,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "LM"))
# Carregar o pacote broom
library(broom)
tidy(fit_lm) |> arrange(p.value) |> filter(p.value < 0.05)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 age 1.66 0.179 9.27 1.94e-20
## 2 flag_X0 -1.45 0.234 -6.19 6.14e-10
## 3 avg_glucose_level 0.200 0.0666 3.01 2.60e- 3
## 4 hypertension_X0 -0.414 0.186 -2.23 2.55e- 2Criação e ajuste dos modelos
Elastic Net - sem interação
Interação das alterações dos hiperparametros lambda e alpha na AUC
# Definir os hiperparâmetros para tunagem
enet <- logistic_reg(penalty = tune(), mixture = tune()) %>%
set_engine("glmnet") %>%
set_mode("classification")
# CV
set.seed(234)
cv_split <- vfold_cv(treinamento, v = 10)
# Paralelização
doParallel::registerDoParallel(makeCluster(16))
tempo <- system.time({
enet_grid <- tune_grid(enet,
receita,
resamples = cv_split,
grid = 100,
metrics = metric_set(roc_auc))
})
autoplot(enet_grid) # Plota os resultados do grid search
Melhor grid de parâmetros de acordo com a métrica AUC
enet_grid %>%
collect_metrics() %>% # Visualiza o tibble de resultados
arrange(desc(mean)) |> head()## # A tibble: 6 × 8
## penalty mixture .metric .estimator mean n std_err .config
## <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.0146 0.482 roc_auc binary 0.862 10 0.0106 Preprocessor1_Model046
## 2 0.00491 0.631 roc_auc binary 0.862 10 0.0114 Preprocessor1_Model062
## 3 0.0112 0.223 roc_auc binary 0.862 10 0.0111 Preprocessor1_Model019
## 4 0.00328 0.421 roc_auc binary 0.861 10 0.0114 Preprocessor1_Model040
## 5 0.00580 0.882 roc_auc binary 0.861 10 0.0112 Preprocessor1_Model088
## 6 0.00850 0.779 roc_auc binary 0.861 10 0.0109 Preprocessor1_Model077Variáveis mais importantes pela biblioteca VIP (que mais afetam o erro do modelo caso sejam alteradas aleatóriamente)
# Seleciona o melhor conjunto de hiperparâmetros
best_enet <- enet_grid %>%
select_by_one_std_err("roc_auc") #select_by_one_std_err() ou select_best()
# tempo #grid = 50, 106s
# Faz o fit do modelo com o melhor conjunto de hiperparâmetros
enet_fit <- finalize_model(enet, parameters = best_enet) %>%
fit(stroke ~ ., data = treinamento_proc)
pred_enet <- predict(enet_fit, new_data = teste_proc, type = "prob")$.pred_1 #pega coluna que dá probabilidade de ter o derrame
desempenho <- desempenho %>%
bind_rows(tibble(prob = pred_enet,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "Elastic Net"))
#desempenho %>% view()
# visualizando variáveis mais importantes
vip(enet_fit, num_features = 20)
Variáveis mais importantes: de maior beta em módulo
# Carregar o pacote broom
library(broom)
# Obter os coeficientes e o intercepto em ordem decrescente
tidy(enet_fit) |> arrange(desc(abs(estimate))) |> head(8)## # A tibble: 8 × 3
## term estimate penalty
## <chr> <dbl> <dbl>
## 1 (Intercept) -3.36 0.00000207
## 2 age 1.62 0.00000207
## 3 work_type_children 1.58 0.00000207
## 4 work_type_outros -1.36 0.00000207
## 5 flag_X0 -0.728 0.00000207
## 6 flag_outros 0.712 0.00000207
## 7 work_type_Self.employed -0.457 0.00000207
## 8 cluster_X2 -0.395 0.00000207Modelo Random Forest
Tunning dos hiperparâmetros
# Definir os hiperparâmetros para tunagem
rf <- rand_forest(trees = tune(), mtry = tune(), min_n = tune()) %>%
set_engine("ranger", importance = "impurity") %>%
set_mode("classification")
# CV
set.seed(234)
cv_split <- vfold_cv(treinamento, v = 5)
# Paralelização
doParallel::registerDoParallel(makeCluster(16))
tempo <- system.time({
rf_grid <- tune_grid(rf,
receita,
resamples = cv_split,
grid = 50,
metrics = metric_set(roc_auc))
})
autoplot(rf_grid) # Plota os resultados do grid search
Melhor grid de parâmetros de acordo com a métrica AUC
rf_grid %>%
collect_metrics() %>% # Visualiza o tibble de resultados
arrange(desc(mean))## # A tibble: 50 × 9
## mtry trees min_n .metric .estimator mean n std_err .config
## <int> <int> <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 3 1061 34 roc_auc binary 0.850 5 0.0110 Preprocessor1_Model…
## 2 1 1313 24 roc_auc binary 0.850 5 0.0124 Preprocessor1_Model…
## 3 3 1018 37 roc_auc binary 0.849 5 0.0116 Preprocessor1_Model…
## 4 2 1592 21 roc_auc binary 0.849 5 0.0110 Preprocessor1_Model…
## 5 2 1486 6 roc_auc binary 0.848 5 0.0109 Preprocessor1_Model…
## 6 8 1414 37 roc_auc binary 0.846 5 0.0109 Preprocessor1_Model…
## 7 10 536 35 roc_auc binary 0.846 5 0.0111 Preprocessor1_Model…
## 8 8 1879 31 roc_auc binary 0.845 5 0.0109 Preprocessor1_Model…
## 9 4 1793 25 roc_auc binary 0.845 5 0.0101 Preprocessor1_Model…
## 10 15 570 39 roc_auc binary 0.845 5 0.0108 Preprocessor1_Model…
## # ℹ 40 more rowsVariáveis mais importantes pela biblioteca VIP (que mais afetam o erro do modelo caso sejam alteradas aleatóriamente)
# Seleciona o melhor conjunto de hiperparâmetros
best_rf <- rf_grid %>%
select_best("roc_auc")
#tempo #grid = 50, 106s
# Faz o fit do modelo com o melhor conjunto de hiperparâmetros
rf_fit <- finalize_model(rf, parameters = best_rf) %>%
fit(stroke ~ ., data = treinamento_proc)
pred_rf <- predict(rf_fit, new_data = teste_proc, type = "prob")$.pred_1 #pega coluna que dá probabilidade de ter o derrame
desempenho <- desempenho %>%
bind_rows(tibble(prob = pred_rf,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "Random Forest"))
#desempenho %>% view()
# visualizando variáveis mais importantes
vip(rf_fit, aesthetics = list(fill = "#FF5757")) #melhor para explicar causalidade
Modelo XGBOOST
Tunning dos hiperparâmetros
boost <- boost_tree(trees = tune(), learn_rate = tune(), mtry = tune(),
tree_depth = tune(), min_n = tune(), sample_size = tune()) %>%
set_engine("xgboost") %>%
set_mode("classification")
# cv
set.seed(234)
cv_split <- vfold_cv(treinamento, v = 5)
# otimização de hiperparametro
doParallel::registerDoParallel(16) #colocar nro de cores da maquina
tempo <- system.time({
boost_grid <- tune_grid(boost,
receita,
resamples = cv_split,
grid = 50,
metrics = metric_set(roc_auc))
})
autoplot(boost_grid) # plota os resultados do grid search
Melhor grid de parâmetros de acordo com a métrica AUC
boost_grid %>%
collect_metrics() %>% # visualiza o tibble de resultados
arrange(desc(mean))## # A tibble: 50 × 12
## mtry trees min_n tree_depth learn_rate sample_size .metric .estimator mean
## <int> <int> <int> <int> <dbl> <dbl> <chr> <chr> <dbl>
## 1 5 67 6 9 0.0192 0.911 roc_auc binary 0.864
## 2 24 570 7 2 0.00670 0.739 roc_auc binary 0.862
## 3 3 513 13 14 0.00708 0.946 roc_auc binary 0.861
## 4 15 1115 2 3 0.00419 0.424 roc_auc binary 0.861
## 5 12 1038 21 4 0.00376 0.948 roc_auc binary 0.860
## 6 19 320 19 11 0.0210 0.720 roc_auc binary 0.859
## 7 25 743 4 2 0.00960 0.332 roc_auc binary 0.859
## 8 10 1529 23 6 0.00131 0.834 roc_auc binary 0.858
## 9 23 1876 10 13 0.00106 0.363 roc_auc binary 0.858
## 10 21 1934 20 9 0.00617 0.649 roc_auc binary 0.855
## # ℹ 40 more rows
## # ℹ 3 more variables: n <int>, std_err <dbl>, .config <chr>Variáveis mais importantes pela biblioteca VIP (que mais afetam o erro do modelo caso sejam alteradas aleatóriamente)
(best_xgb <- boost_grid %>%
select_best("roc_auc")) # salva o melhor conjunto de parametros## # A tibble: 1 × 7
## mtry trees min_n tree_depth learn_rate sample_size .config
## <int> <int> <int> <int> <dbl> <dbl> <chr>
## 1 5 67 6 9 0.0192 0.911 Preprocessor1_Model08# tempo #grid = 50, 161s
#faz o fit do modelo com o melhor conjunto de hiperparâmetros
boost_fit <- finalize_model(boost, parameters = best_xgb) %>%
fit(stroke ~ ., treinamento_proc)
pred_xgbm <- predict(boost_fit, new_data = teste_proc, type = "prob")$.pred_1
desempenho <- desempenho %>%
bind_rows(tibble(prob = pred_xgbm,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "XGBM"))
# visualizando variáveis mais importantes
vip(boost_fit)
Redes Neurais
Definição de estrutura e loading de pesos treinados
#desempenho <- desempenho |> filter(metodo != "Neural Networks")
X_trn <- treinamento_proc %>% dplyr::select(-stroke) %>% as.matrix()
X_tst <- teste_proc %>% dplyr::select(-stroke) %>% as.matrix()
y_trn <- treinamento_proc$stroke %>% to_categorical()
y_tst <- teste_proc$stroke %>% to_categorical()
#net <- NULL
# Modelo Sequencial
net <- keras_model_sequential() %>%
layer_dense(units = 512, activation = "relu", input_shape = ncol(X_trn)) %>%
layer_dropout(rate = 0.4) %>%
layer_dense(units = 256, activation = "relu") %>%
layer_dropout(rate = 0.4) %>%
layer_dense(units = 128, activation = "relu") %>%
layer_dropout(rate = 0.3) %>%
layer_dense(units = 64, activation = "relu") %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 32, activation = "relu") %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 16, activation = "relu") %>%
layer_dropout(rate = 0.1) %>%
layer_dense(units = 2, activation = "softmax") # Camada de saída com ativação sigmoid
# definicao do estimador da rede neural
net %>%
compile(loss = "categorical_crossentropy",
optimizer = optimizer_rmsprop(),
metrics = c("Recall")
)
#accuracy Recall(sensibilidade) Precision
summary(net)## Model: "sequential"
## ________________________________________________________________________________
## Layer (type) Output Shape Param #
## ================================================================================
## dense_6 (Dense) (None, 512) 13824
## dropout_5 (Dropout) (None, 512) 0
## dense_5 (Dense) (None, 256) 131328
## dropout_4 (Dropout) (None, 256) 0
## dense_4 (Dense) (None, 128) 32896
## dropout_3 (Dropout) (None, 128) 0
## dense_3 (Dense) (None, 64) 8256
## dropout_2 (Dropout) (None, 64) 0
## dense_2 (Dense) (None, 32) 2080
## dropout_1 (Dropout) (None, 32) 0
## dense_1 (Dense) (None, 16) 528
## dropout (Dropout) (None, 16) 0
## dense (Dense) (None, 2) 34
## ================================================================================
## Total params: 188946 (738.07 KB)
## Trainable params: 188946 (738.07 KB)
## Non-trainable params: 0 (0.00 Byte)
## ________________________________________________________________________________Treina o modelo de redes neurais e salva
tensorflow::set_random_seed(123)
# Definindo o callback de Early Stopping
callback <- callback_early_stopping(
monitor = "val_loss", # Métrica a ser monitorada
patience = 50, # Número de épocas para esperar após a parada de melhorias
restore_best_weights = TRUE) # Restaura os pesos do modelo para a melhor época
history <- net %>%
fit(X_trn, y_trn, epochs = 10000,
batch_size = 80, validation_split = 0.2, #DIVIDE O TOTAL DE DADOS EM 80 BLOCOS
callbacks = list(callback)) #pega 20% pra validar enquanto roda
#plot(history)
# Salvando os pesos do modelo
save_model_weights_hdf5(net, "pesos_nn_1.h5")carrega pesos salvos, prediz e carrega predição na base comparativa
# Carregando os pesos salvos no modelo
net %>% load_model_weights_hdf5("pesos_nn_1.h5")
pred_nn <- predict(net, X_tst)[,2] #prob de apresentar o evento## 32/32 - 0s - 121ms/epoch - 4ms/stepdesempenho <- desempenho %>%
bind_rows(tibble(prob = pred_nn,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "Neural Networks"))Modelo de redes neurais com 200 epocas (sem call back)
Treina o modelo de redes neurais e salva
tensorflow::set_random_seed(42)
history <- net %>%
fit(X_trn, y_trn, epochs = 200,
batch_size = 80, validation_split = 0.2)
# Salvando os pesos do modelo
save_model_weights_hdf5(net, "pesos_nn_2.h5")# Carregando os pesos salvos no modelo
net %>% load_model_weights_hdf5("pesos_nn_2.h5")
pred_nn_2 <- predict(net, X_tst)[,2] #prob de apresentar o evento## 32/32 - 0s - 84ms/epoch - 3ms/stepdesempenho <- desempenho %>%
bind_rows(tibble(prob = pred_nn_2,
classes = ifelse(teste_proc$stroke == "1", "Yes", "No"),
metodo = "Neural Networks - sem callback"))Comparativo de todos os modelos
Como esperado redes neurais com callback teve melhor desempenho que sem callback
desempenho %>%
mutate(classes = factor(classes)) %>%
group_by(metodo) %>%
roc_auc(classes, prob, event_level = "second") %>%
arrange(desc(.estimate))## # A tibble: 8 × 4
## metodo .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 Elastic Net roc_auc binary 0.823
## 2 LM roc_auc binary 0.822
## 3 LASSO - interact roc_auc binary 0.820
## 4 Neural Networks roc_auc binary 0.818
## 5 XGBM roc_auc binary 0.803
## 6 Random Forest roc_auc binary 0.799
## 7 LM - interact roc_auc binary 0.782
## 8 Neural Networks - sem callback roc_auc binary 0.767Curva AUC
desempenho %>%
mutate(classes = factor(classes)) %>%
group_by(metodo) %>%
roc_curve(classes, prob, event_level = "second") %>%
autoplot()
Análise aprofundada dos modelos e pontos de corte
Para mensurar o melhor ponto de corte definimos uma variavel chamada “campo” que é o produto de ppv por sensibilidade. Para dar uma importancia maior para sensibilidade elevamos ela ao cubo. Para ficar mais fácil de vizualisar o resultado final multiplicamos por 100. Ordenamos em ordem decrescente os 3 principais pontos de corte. Estamos considerando todos os gráficos pois, não necessariamente o modelo com maior auc será o melhor para nosso problema em questão.
accuracy : acertos/casos totais
sensitivity/recall : acertos positivos/total de positivos
specificity : assertos negativos/total de negativos reais 1-specificity : erros negativos/total de negativos reais (% dos totais que não terão derrame que fizeram exame em vão)
ppv/precisão : assertos positivos/total de positivos previstos (% das pessoas chamadas que realmente terão derrame) npv : assertos negativos/total de negativos previstos (% das pessoas que não serão chamadas que realmente não terão derrame)
recall : assertos positivos/total de positivos reais
Crio função para calcular formula de campo para cada modelo para cada ponto de corte utilizando calculo que leve em consideração os pesos da sensibilidade e especificidade para o negócio em questão. Melhor modelo e melhor ponto de corte obtido:
# Função para calcular o valor do campo
calc_campo <- function(pred, model_name) {
coords(roc(teste_proc$stroke, pred), c(seq(.001,.5,.001)),
ret = c("threshold","accuracy", "sensitivity", "specificity", "1-specificity", "ppv", "npv")) %>%
mutate(campo = (sensitivity^3)*specificity, model = model_name) %>%
arrange(desc(campo))
}
# Lista de previsões
pred_list <- list(LASSO = pred_LASSO, lm = pred_lm, enet = pred_enet, rf = pred_rf, xgbm = pred_xgbm, nn = pred_nn, nnscall = pred_nn_2)
# Aplica a função a cada conjunto de previsões
campo_values <- lapply(names(pred_list), function(x) calc_campo(pred_list[[x]], x))
# Combina os resultados em um data frame
campo_df <- do.call(rbind, campo_values)
# Retorna a linha com o maior valor do campo
campo_df[which.max(campo_df$campo), ]## threshold accuracy sensitivity specificity 1-specificity ppv npv campo
## 1 0.021 0.5372 0.963 0.5134 0.4866 0.09943 0.996 0.4585
## model
## 1 LASSOMatriz de confusão do melhor modelo e ponto de corte
corte <- 0.021
previsoes_1_1 <- ifelse(pred_LASSO >= corte, 1, 0) %>% as.tibble
previsoes_1_1 <- previsoes_1_1 %>% bind_cols(tibble(obs = teste_proc$stroke))
table(previsoes_1_1) |> `dimnames<-`(list(previsto = c("0", "1"), real = c("0", "1")))## real
## previsto 0 1
## 0 497 2
## 1 471 52Conclusão: ao fazer exame nos casos com 2.1% ou mais de probabilidade segundo previsão do modelo LASSO, estaremos diminuindo em aproximadamente 46% a quantidade de exames necessárias enquanto somente 2 pessoas com derrame ficariam de fora da previsão. Isso gerara ganhos consideraveis para o hospital e ganhos indiretos aos pacientes, uma vez que com menos pacientes na fila os exames teoricamente devem sair mais rápido e o atendimento deve ser melhor.
De que forma podemos obter ganhos diretos para os pacientes? E se quisermos dar prioridade para casos com probabilidade de derrame mais elevada?
pro ponto de corte de 20% qual modelo possui melhor valor de “campo” ?
xgbm é o melhor modelo para um ponto de corte de 20% de probabilidade de derrame
#escolhe a resposta com maior campo. desta vez sensibilidade (quem é deixado de fora) tem peso normal pois step anterior já leva ela em consideração
# Função para calcular o valor do campo
calc_campo <- function(pred) {
coords(roc(teste_proc$stroke, pred), .2,
ret = c("threshold","accuracy", "sensitivity", "specificity", "1-specificity", "ppv", "npv")) %>%
mutate(campo = (sensitivity^1)*specificity*10) %>%
pull(campo)
}
# Lista de previsões
pred_list <- list(LASSO = pred_LASSO, lm = pred_lm, enet = pred_enet, rf = pred_rf, xgbm = pred_xgbm, nn = pred_nn, nnscall = pred_nn_2)
# Aplica a função a cada conjunto de previsões
(campo_values <- lapply(pred_list, calc_campo))## $LASSO
## [1] 1.081
##
## $lm
## [1] 1.61
##
## $enet
## [1] 1.61
##
## $rf
## [1] 2.169
##
## $xgbm
## [1] 4.637
##
## $nn
## [1] 1.422
##
## $nnscall
## [1] 2.228Lista de prioridade de exames obtida após executar os modelos:
corte <- 0.2 #20%
previsoes_2 <- ifelse(pred_xgbm >= corte, 1, 0) %>% as.tibble %>% setNames("valor2")
lista_prioritaria <- teste_proc %>%
bind_cols(previsoes_1_1) %>% bind_cols(previsoes_2)
lista_prioritaria <- lista_prioritaria %>% mutate(resultado=ifelse(value == 1, "em risco", "sem risco"),
resultado = ifelse(valor2 == 1, "urgente", resultado)) %>% dplyr::select(-value, -valor2)
lista_prioritaria$resultado %>% table()## .
## em risco sem risco urgente
## 333 499 190Percentual de pacientes que teriam derrame e ficariam de fora (falsos negativos):
(lista_prioritaria %>% filter(resultado == 'sem risco' & stroke == 1) %>% nrow()/lista_prioritaria %>% dplyr::filter(stroke == 1) %>% nrow())*100## [1] 3.704Percentual de exames diminuidos:
#Percentual de exames diminuidos
(1-(lista_prioritaria %>% filter(resultado != 'sem risco') %>% nrow()/nrow(previsoes_1_1)))*100## [1] 48.83Geramos uma lista indicando se o paciente precisará ou não realizar um exame e, caso necessário, classificamos como “urgente” aqueles com probabilidade igual ou superior a 20% de apresentar um derrame.
Curiosamente, um dos modelos com AUC mais baixa em comparação com os outros (XGBM) foi o mais eficaz na previsão de casos críticos (probabilidade de AVC acima de 20%). Isso demonstra que, para modelos de classificação, é mais importante ter um modelo específico para o problema em questão do que um modelo com maior AUC para todos os pontos de corte.