library(caret)
library(tidyverse)
library(neuralnet)Base de Seguros
um exemplo de RNA para classificação
O Problema
Considere a base de dados insurance.csv, composta pelas variáveis: age, sex, bmi, children, smoker, region, married, charges. Para que possamos trabalhar com essa base em modelos de classificação vamos substituir a variável alvo charge para uma variável categórica que indica se o cliente é de custo elevado ou não.
Vamos ajustar diferentes modelos de classificação para prever a categoria do cliente, custo alto ou custo baixo, a partir das suas características como idade, gênero, índice de massa corporal, número de filhos, hábito de fumar, região de residência e se é casado ou não.
Pacotes necessários
O primeiro passo será carregar os pacotes necessários.
Leitura da base de dados
Aqui vou aproveitar o que já foi feito para o problema de regressão e carregar a base sem NAs antes de passar pela padronização.
base_treino = readRDS("base_treino_sem_NA.RDS")
glimpse(base_treino)Rows: 936
Columns: 7
$ age <dbl> 19, 18, 28, 32, 31, 37, 37, 60, 25, 62, 56, 27, 19, 52, 30, 6…
$ sex <fct> female, male, male, male, female, female, male, female, male,…
$ bmi <dbl> 27.900, 33.770, 33.000, 28.880, 25.740, 27.740, 29.830, 25.84…
$ children <dbl> 0, 1, 3, 0, 0, 3, 2, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 2, 1, 0…
$ smoker <fct> yes, no, no, no, no, no, no, no, no, yes, no, yes, no, no, ye…
$ region <fct> southwest, southeast, southeast, northwest, southeast, northw…
$ charges <dbl> 16884.924, 1725.552, 4449.462, 3866.855, 3756.622, 7281.506, …A criação da nova variável alvo será feita da seguinte maneira:
hist(base_treino$charges)
custo = ifelse(base_treino$charges>20000,"alto","normal")
base_treino = base_treino |> cbind(custo)
base_treino = base_treino |> select(-charges)Padronização das variáveis independentes numéricas
As variáveis independentes (aquelas que não são a variável resposta) numéricas devem ser padronizadas e podem ser padronizadas pela forma \[ \tilde{X} = \dfrac{X - \bar{X}}{sd(X)} \] E esta padronização pode ser realizada com a função scale.
age_ = scale(base_treino$age)
attr(age_,"scaled:center")[1] 39.25321attr(age_,"scaled:scale")[1] 14.19539bmi_ = scale(base_treino$bmi)
children_ = scale(base_treino$children)Tratamento nas variáveis categóricas
As variáveis categóricas precisam ser transformadas em indicadoras antes do ajuste de uma rede neural pela função neuralnet. Para isso vamos primeiro criar uma base de dados com as variáveis numéricas padronizadas e as variáveis categóricas como factor. Depois usamos a função modelmatrix e por último eliminamos a coluna intersept.
base_treino_ = tibble(
age = as.numeric(age_),
bmi = as.numeric(bmi_),
children = as.numeric(children_),
sex = base_treino$sex,
smoker = base_treino$smoker,
region = base_treino$region,
custo = base_treino$custo
)
glimpse(base_treino_)Rows: 936
Columns: 7
$ age <dbl> -1.4267453, -1.4971907, -0.7927366, -0.5109550, -0.5814004, -…
$ bmi <dbl> -0.45420835, 0.49379114, 0.36943687, -0.29593927, -0.80304632…
$ children <dbl> -0.92254085, -0.08337833, 1.59494671, -0.92254085, -0.9225408…
$ sex <fct> female, male, male, male, female, female, male, female, male,…
$ smoker <fct> yes, no, no, no, no, no, no, no, no, yes, no, yes, no, no, ye…
$ region <fct> southwest, southeast, southeast, northwest, southeast, northw…
$ custo <chr> "normal", "normal", "normal", "normal", "normal", "normal", "…base_treino_final =
model.matrix( ~ age + sex + bmi +
children + smoker +
region + custo,
data = base_treino_)
head(base_treino_final) (Intercept) age sexmale bmi children smokeryes
1 1 -1.4267453 0 -0.4542083 -0.92254085 1
2 1 -1.4971907 1 0.4937911 -0.08337833 0
3 1 -0.7927366 1 0.3694369 1.59494671 0
4 1 -0.5109550 1 -0.2959393 -0.92254085 0
5 1 -0.5814004 0 -0.8030463 -0.92254085 0
6 1 -0.1587280 0 -0.4800482 1.59494671 0
regionnorthwest regionsoutheast regionsouthwest custonormal
1 0 0 1 1
2 0 1 0 1
3 0 1 0 1
4 1 0 0 1
5 0 1 0 1
6 1 0 0 1base_treino_final = base_treino_final[,-1]
head(base_treino_final) age sexmale bmi children smokeryes regionnorthwest
1 -1.4267453 0 -0.4542083 -0.92254085 1 0
2 -1.4971907 1 0.4937911 -0.08337833 0 0
3 -0.7927366 1 0.3694369 1.59494671 0 0
4 -0.5109550 1 -0.2959393 -0.92254085 0 1
5 -0.5814004 0 -0.8030463 -0.92254085 0 0
6 -0.1587280 0 -0.4800482 1.59494671 0 1
regionsoutheast regionsouthwest custonormal
1 0 1 1
2 1 0 1
3 1 0 1
4 0 0 1
5 1 0 1
6 0 0 1Ajuste do Modelo de Rede Neural com um único neurônio
A função usada para isso será a neuralnet do pacote neuralnet.
Modelo Completo função de ativação logística
modelo_completo_log = neuralnet(
formula = custonormal ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast + regionsouthwest,
data = base_treino_final,
hidden = 0,
linear.output = FALSE,
err.fct = "ce")Podemos visualizar o modelo criado com a função plot:
plot(modelo_completo_log,rep="best")
Podemos acessar os valores estimados pelos pesos sinápticos da seguinte forma:
modelo_completo_log$weights[[1]][[1]][,1][1] 2.92611918 -0.69307925 -0.17821346 -0.53749487 0.03676583 -4.44715032
[7] 0.09914370 0.19944378 0.35317416Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:
head(modelo_completo_log$net.result[[1]]) [,1]
1 0.5078012
2 0.9762595
3 0.9663284
4 0.9653211
5 0.9806635
6 0.9692889Faça na mão a estimativa para a primeira observação da base e verifique se está correta:
pesos = modelo_completo_log$weights[[1]][[1]][,1]
x = as.numeric(modelo_completo_log$data[1,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest")])
x = c(1,x)
soma = sum(x*pesos)
1/(1+exp(-soma))[1] 0.5078012modelo_completo_log$net.result[[1]][1][1] 0.5078012Modelo sem Região com função de ativação logística
modelo_sem_regiao_log = neuralnet(
formula = custonormal ~ age + sexmale + bmi + children + smokeryes,
data = base_treino_final,
hidden = 0,
linear.output = FALSE,
err.fct = "ce")Podemos visualizar o modelo criado com a função plot:
plot(modelo_sem_regiao_log,rep="best")
Podemos acessar os valores estimados pelos pesos sinápticos da seguinte forma:
modelo_sem_regiao_log$weights[[1]][[1]][,1][1] 3.07758608 -0.69634809 -0.16472247 -0.51571932 0.04460086 -4.43638983Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:
head(modelo_sem_regiao_log$net.result[[1]]) [,1]
1 0.4570574
2 0.9758035
3 0.9659570
4 0.9670781
5 0.9792744
6 0.9708813Faça na mão a estimativa para a primeira observação da base e verifique se está correta:
pesos = modelo_sem_regiao_log$weights[[1]][[1]][,1]
x = as.numeric(modelo_sem_regiao_log$data[1,c("age","sexmale","bmi","children","smokeryes")])
x = c(1,x)
soma = sum(x*pesos)
(1)/(1+exp(-soma))[1] 0.4570574modelo_sem_regiao_log$net.result[[1]][1][1] 0.4570574Medidas de Qualidade na base de treino - problemas de Classificação
Vamos comparar o desempenho dos dois modelos de RNA ajustados para prever o a categoria do curso do cliente. Para isso, primeiro, vamos guardar em um objeto do R os valores previstos.
y_completo_log = modelo_completo_log$net.result[[1]]
y_sem_regiao_log = modelo_sem_regiao_log$net.result[[1]]Entropia Cruzada
A Entropia Cruzada é uma medida de desempenho de modelos de classificação que avalia o quanto bom está o valor previsto \(\hat{y}_i\), antes mesmo do pesquisador ter definido a classe prevista.
\[ \begin{array}{ll} EC &= - \sum_{i=1}^N \left( y_i\ln(\hat{y}_i) + (1-y_i)\ln(1-\hat{y}_i) \right) \\ &= - \sum_{i=1}^N EC_i \quad \text{, sendo } EC_i = y_i\ln(\hat{y}_i) + (1-y_i)\ln(1-\hat{y}_i). \end{array} \]
No R podemos calculá-la assim:
y = base_treino_final[,"custonormal"]
(EC_completo_log = -sum(y*log(y_completo_log) + (1-y)*log(1-y_completo_log)))[1] 255.2122(EC_sem_regiao_log = -sum(y*log(y_sem_regiao_log) + (1-y)*log(1-y_sem_regiao_log)))[1] 255.7868Curva ROC e Matriz de Confusão
Primeiro vamos analisar a curva ROC para o modelo completo e definir o valor de corte para realizar as previsões das classes.
library(pROC)
roc_completo = roc(response=y,
predictor=y_completo_log)
plot(roc_completo)
auc(roc_completo)Area under the curve: 0.8827A partir do obejto retornado pela função roc podemos encontrar o valor ótimo de corte para definir o que será classificado como 1 (normal) e como 0 (alto).
(otimo = coords(roc_completo, "best", best.method = "youden")) threshold specificity sensitivity
1 0.4036661 0.7305699 0.9771198c_completo = otimo$threshold
prev_completo = ifelse(y_completo_log < c_completo,0,1)Agora podemos usar a função confusionMatrix do pacote caret para definir a matriz de confusão e as medidas de qualidade geradas a partir dela.
CM_completo = confusionMatrix(
factor(prev_completo),
factor(y)
)
CM_completoConfusion Matrix and Statistics
Reference
Prediction 0 1
0 141 17
1 52 726
Accuracy : 0.9263
95% CI : (0.9076, 0.9422)
No Information Rate : 0.7938
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.7586
Mcnemar's Test P-Value : 4.256e-05
Sensitivity : 0.7306
Specificity : 0.9771
Pos Pred Value : 0.8924
Neg Pred Value : 0.9332
Prevalence : 0.2062
Detection Rate : 0.1506
Detection Prevalence : 0.1688
Balanced Accuracy : 0.8538
'Positive' Class : 0
Faremos o mesmo para o modelo sem região.
roc_sem_regiao = roc(response=y,
predictor=y_sem_regiao_log)
plot(roc_sem_regiao)
auc(roc_sem_regiao)Area under the curve: 0.8821(otimo = coords(roc_sem_regiao, "best", best.method = "youden")) threshold specificity sensitivity
1 0.3892133 0.7253886 0.9811575c_sem_regiao = otimo$threshold
prev_sem_regiao = ifelse(y_sem_regiao_log < c_sem_regiao,0,1)CM_sem_regiao = confusionMatrix(
factor(prev_sem_regiao),
factor(y)
)
CM_sem_regiaoConfusion Matrix and Statistics
Reference
Prediction 0 1
0 140 14
1 53 729
Accuracy : 0.9284
95% CI : (0.91, 0.9441)
No Information Rate : 0.7938
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.7637
Mcnemar's Test P-Value : 3.443e-06
Sensitivity : 0.7254
Specificity : 0.9812
Pos Pred Value : 0.9091
Neg Pred Value : 0.9322
Prevalence : 0.2062
Detection Rate : 0.1496
Detection Prevalence : 0.1645
Balanced Accuracy : 0.8533
'Positive' Class : 0
Medidas de Qualidade na base de teste - problemas de Classificação
base_teste = readRDS("base_teste_bruta.RDS")
glimpse(base_teste)Rows: 400
Columns: 8
$ age <dbl> 33, 46, 23, 23, 56, 30, 37, 59, 63, 31, 22, 18, 19, 63, 28, 6…
$ sex <fct> male, female, male, male, male, female, male, female, female,…
$ bmi <dbl> 22.705, 33.440, 34.400, 23.845, 40.300, 32.400, 28.025, 27.72…
$ children <dbl> 0, 1, 0, 0, 0, 1, 2, 3, 0, 2, 0, 0, 5, 0, 1, 3, 1, 0, 1, 3, 2…
$ smoker <fct> no, no, no, no, no, no, no, no, no, yes, yes, no, no, no, yes…
$ region <fct> northwest, southeast, southwest, northeast, southwest, southw…
$ charges <dbl> 21984.471, 8240.590, 1826.843, 2395.172, 10602.385, NA, 6203.…
$ married <fct> no, yes, NA, yes, yes, NA, NA, NA, yes, yes, yes, no, NA, yes…Transformações na Base de Teste
Para calcular as medidas de qualidades do ajuste na base de teste, primeiro é preciso realizar na base de teste as mesmas transformações feitas na base de treino, que para esse exemplo foram: retiradas das linhas com NAs e transforações nas variáveis.
base_teste = base_teste |> select(-married)
base_teste = na.omit(base_teste)
m_age = attr(age_,"scaled:center")
s_age = attr(age_,"scaled:scale")
age_teste = (base_teste$age - m_age)/s_age
m_bmi = attr(bmi_,"scaled:center")
s_bmi = attr(bmi_,"scaled:scale")
bmi_teste = (base_teste$bmi - m_bmi)/s_bmi
m_children = attr(children_,"scaled:center")
s_children = attr(children_,"scaled:scale")
children_teste = (base_teste$children - m_children)/s_childrenAgora vamos criar a base de teste transformadae realizar o tratamento nas variáveis categóricas.
base_teste_ = tibble(
age = age_teste,
bmi = bmi_teste,
children = children_teste,
sex = base_teste$sex,
smoker = base_teste$smoker,
region = base_teste$region
)
base_teste_final = model.matrix( ~ age + sex + bmi + children + smoker + region,
data = base_teste_)
base_teste_final = base_teste_final[,-1]É importante que as variáveis categóricas criadas pela função model.matrix sejam as mesmas que foram criadas para a base de treino. Vamos verificar.
colnames(base_treino_final)[1] "age" "sexmale" "bmi" "children"
[5] "smokeryes" "regionnorthwest" "regionsoutheast" "regionsouthwest"
[9] "custonormal" colnames(base_teste_final)[1] "age" "sexmale" "bmi" "children"
[5] "smokeryes" "regionnorthwest" "regionsoutheast" "regionsouthwest"Previsão do modelo
Agora a base de teste está pronta para realizarmos as previsões, que virão na unidade transformada.
pred_completo_log = predict(modelo_completo_log,newdata = base_teste_final)[,1]
pred_sem_regiao_log = predict(modelo_sem_regiao_log,newdata = base_teste_final)[,1]Medidas de Qualidade
Primeiro vamos guardar o valor real da variável alvo na base de teste.
y = ifelse(base_teste$charges > 20000,0,1)O valor de corte utilizado faz parte do modelo, logo ele será aquele definido no treino.
prev_completo_teste = ifelse(pred_completo_log < c_completo,0,1)
prev_sem_regiao_teste = ifelse(pred_sem_regiao_log < c_sem_regiao,0,1)CM_completo_teste = confusionMatrix(
factor(prev_completo_teste),
factor(y)
)
CM_completo_testeConfusion Matrix and Statistics
Reference
Prediction 0 1
0 66 9
1 12 310
Accuracy : 0.9471
95% CI : (0.9203, 0.967)
No Information Rate : 0.8035
P-Value [Acc > NIR] : <2e-16
Kappa : 0.83
Mcnemar's Test P-Value : 0.6625
Sensitivity : 0.8462
Specificity : 0.9718
Pos Pred Value : 0.8800
Neg Pred Value : 0.9627
Prevalence : 0.1965
Detection Rate : 0.1662
Detection Prevalence : 0.1889
Balanced Accuracy : 0.9090
'Positive' Class : 0
CM_sem_regiao_teste = confusionMatrix(
factor(prev_sem_regiao_teste),
factor(y)
)
CM_sem_regiao_testeConfusion Matrix and Statistics
Reference
Prediction 0 1
0 66 12
1 12 307
Accuracy : 0.9395
95% CI : (0.9114, 0.9609)
No Information Rate : 0.8035
P-Value [Acc > NIR] : 1.01e-14
Kappa : 0.8085
Mcnemar's Test P-Value : 1
Sensitivity : 0.8462
Specificity : 0.9624
Pos Pred Value : 0.8462
Neg Pred Value : 0.9624
Prevalence : 0.1965
Detection Rate : 0.1662
Detection Prevalence : 0.1965
Balanced Accuracy : 0.9043
'Positive' Class : 0