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
Em seguida faremos a leitura da base de dados. Nesse momento é preciso verificar o tipo das variáveis e todas as variáveis categóricas precisam ser transformadas em factor.
base = read_csv("insurance.csv")
glimpse(base)Rows: 1,338
Columns: 9
$ age <dbl> 19, 18, 28, 33, 32, 31, 46, 37, 37, 60, 25, 62, 23, 56, 27, 1…
$ sex <chr> "female", "male", "male", "male", "male", "female", "female",…
$ bmi <dbl> 27.900, 33.770, 33.000, 22.705, 28.880, 25.740, 33.440, 27.74…
$ children <dbl> 0, 1, 3, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0…
$ smoker <chr> "yes", "no", "no", "no", "no", "no", "no", "no", "no", "no", …
$ region <chr> "southwest", "southeast", "southeast", "northwest", "northwes…
$ charges <dbl> 16884.924, 1725.552, 4449.462, 21984.471, 3866.855, 3756.622,…
$ married <chr> "yes", NA, "yes", "no", "yes", "yes", "yes", "yes", "yes", NA…
$ country <chr> "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa"…base$sex = factor(base$sex)
base$smoker = factor(base$smoker)
base$region = factor(base$region)
base$married = factor(base$married)
base$country = factor(base$country)
glimpse(base)Rows: 1,338
Columns: 9
$ age <dbl> 19, 18, 28, 33, 32, 31, 46, 37, 37, 60, 25, 62, 23, 56, 27, 1…
$ sex <fct> female, male, male, male, male, female, female, female, male,…
$ bmi <dbl> 27.900, 33.770, 33.000, 22.705, 28.880, 25.740, 33.440, 27.74…
$ children <dbl> 0, 1, 3, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0…
$ smoker <fct> yes, no, no, no, no, no, no, no, no, no, no, yes, no, no, yes…
$ region <fct> southwest, southeast, southeast, northwest, northwest, southe…
$ charges <dbl> 16884.924, 1725.552, 4449.462, 21984.471, 3866.855, 3756.622,…
$ married <fct> yes, NA, yes, no, yes, yes, yes, yes, yes, NA, NA, yes, NA, N…
$ country <fct> usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, u…E a criação da nova variável alvo:
hist(base$charges)
custo = ifelse(base$charges>20000,"alto","normal")
base = base |> cbind(custo)Divisão da base em treino e teste
Em seguida dividimos a base em treino e teste de forma que 70% ficará no treino e 30% no teste. Para garantir a reprodutibilidade do código devemos fixar a semente antes de realizar a partição.
set.seed(123456789)
N = dim(base)[1]
indices_treino = createDataPartition(1:N,p=0.7)[[1]]
base_treino = base[indices_treino,]
base_teste = base[-indices_treino,]Tratamento das NAs
A partir de agora vamos usar somente a base de treino. A base de teste só volta ao uso quando quisermos medir a qualidade das previsões fora da amostra.
Sobre o tratamento das NAs, temos a alternativa de eliminar as linhas, de imputar valores ou eliminar colunas. Neste exemplo, vamos simplesmente eliminar as linhas que contém algum NA de forma a simplificar esta etapa, mas nem sempre está é a alternativa mais interessante.
sum(is.na(base_treino)) #num total de NAs na base[1] 384base_treino = na.omit(base_treino)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] 40.09009attr(age_,"scaled:scale")[1] 14.32938bmi_ = 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,
married = base_treino$married,
custo = base_treino$custo
)
glimpse(base_treino_)Rows: 555
Columns: 8
$ age <dbl> -1.471807263, -0.843727188, -0.564580488, -0.634367163, -0.21…
$ bmi <dbl> -0.4625470253, 0.3568875837, -0.3050870416, -0.8096016832, -0…
$ children <dbl> -0.9223682, 1.4999394, -0.9223682, -0.9223682, 1.4999394, 0.6…
$ sex <fct> female, male, male, female, female, male, female, female, mal…
$ smoker <fct> yes, no, no, no, no, no, yes, no, yes, no, no, yes, no, no, y…
$ region <fct> southwest, southeast, northwest, southeast, northwest, northe…
$ married <fct> yes, yes, yes, yes, yes, yes, yes, no, yes, yes, no, yes, yes…
$ custo <chr> "normal", "normal", "normal", "normal", "normal", "normal", "…base_treino_final =
model.matrix( ~ age + sex + bmi +
children + smoker +
region + married + custo,
data = base_treino_)
head(base_treino_final) (Intercept) age sexmale bmi children smokeryes
1 1 -1.4718073 0 -0.4625470 -0.9223682 1
2 1 -0.8437272 1 0.3568876 1.4999394 0
3 1 -0.5645805 1 -0.3050870 -0.9223682 0
4 1 -0.6343672 0 -0.8096017 -0.9223682 0
5 1 -0.2156471 0 -0.4882548 1.4999394 0
6 1 -0.2156471 1 -0.1524473 0.6925035 0
regionnorthwest regionsoutheast regionsouthwest marriedyes custonormal
1 0 0 1 1 1
2 0 1 0 1 1
3 1 0 0 1 1
4 0 1 0 1 1
5 1 0 0 1 1
6 0 0 0 1 1base_treino_final = base_treino_final[,-1]
head(base_treino_final) age sexmale bmi children smokeryes regionnorthwest
1 -1.4718073 0 -0.4625470 -0.9223682 1 0
2 -0.8437272 1 0.3568876 1.4999394 0 0
3 -0.5645805 1 -0.3050870 -0.9223682 0 1
4 -0.6343672 0 -0.8096017 -0.9223682 0 0
5 -0.2156471 0 -0.4882548 1.4999394 0 1
6 -0.2156471 1 -0.1524473 0.6925035 0 0
regionsoutheast regionsouthwest marriedyes custonormal
1 0 1 1 1
2 1 0 1 1
3 0 0 1 1
4 1 0 1 1
5 0 0 1 1
6 0 0 1 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 + marriedyes,
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)Podemos acessar os valores estimados pelos pesos sinápticos da seguinte forma:
modelo_completo_log$weights[[1]][[1]][,1] [1] 2.82248641 -0.72095405 -0.09973425 -0.50826106 0.01322051 -4.50411673
[7] 0.19431828 -0.19291685 0.21281538 0.16413254Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:
head(modelo_completo_log$net.result[[1]]) [,1]
1 0.4948310
2 0.9585405
3 0.9741967
4 0.9746759
5 0.9735167
6 0.9580730Faç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","marriedyes")])
x = c(1,x)
soma = sum(x*pesos)
1/(1+exp(-soma))[1] 0.494831modelo_completo_log$net.result[[1]][1][1] 0.494831Modelo sem Região com função de ativação logística
modelo_sem_regiao_log = neuralnet(
formula = custonormal ~ age + sexmale + bmi + children + smokeryes + marriedyes,
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)Podemos acessar os valores estimados pelos pesos sinápticos da seguinte forma:
modelo_sem_regiao_log$weights[[1]][[1]][,1][1] 2.85705993 -0.72317232 -0.10069231 -0.54147262 0.01246142 -4.52209374
[7] 0.17284912Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:
head(modelo_sem_regiao_log$net.result[[1]]) [,1]
1 0.4529269
2 0.9665869
3 0.9704366
4 0.9804598
5 0.9697906
6 0.9599342Faç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","marriedyes")])
x = c(1,x)
soma = sum(x*pesos)
(1)/(1+exp(-soma))[1] 0.4529269modelo_sem_regiao_log$net.result[[1]][1][1] 0.4529269Medidas 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))
EC_sem_regiao_log = -sum(y*log(y_sem_regiao_log) + (1-y)*log(1-y_sem_regiao_log))
barplot(c(completo_log = EC_completo_log,
sem_regiao_log = EC_completo_log))
Curva 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.8892A 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.4207017 0.7317073 0.974537c_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 90 11
1 33 421
Accuracy : 0.9207
95% CI : (0.895, 0.9418)
No Information Rate : 0.7784
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.7545
Mcnemar's Test P-Value : 0.001546
Sensitivity : 0.7317
Specificity : 0.9745
Pos Pred Value : 0.8911
Neg Pred Value : 0.9273
Prevalence : 0.2216
Detection Rate : 0.1622
Detection Prevalence : 0.1820
Balanced Accuracy : 0.8531
'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.8862(otimo = coords(roc_sem_regiao, "best", best.method = "youden")) threshold specificity sensitivity
1 0.4094996 0.7317073 0.9722222c_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 90 12
1 33 420
Accuracy : 0.9189
95% CI : (0.893, 0.9402)
No Information Rate : 0.7784
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.7497
Mcnemar's Test P-Value : 0.002869
Sensitivity : 0.7317
Specificity : 0.9722
Pos Pred Value : 0.8824
Neg Pred Value : 0.9272
Prevalence : 0.2216
Detection Rate : 0.1622
Detection Prevalence : 0.1838
Balanced Accuracy : 0.8520
'Positive' Class : 0
Medidas de Qualidade na base de teste - problemas de Classificação
glimpse(base_teste)Rows: 400
Columns: 10
$ 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…
$ country <fct> usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, usa, u…
$ custo <chr> "alto", "normal", "normal", "normal", "normal", NA, "normal",…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 = 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,
married = base_teste$married
)
base_teste_final = model.matrix( ~ age + sex + bmi + children + smoker + region + married,
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] "marriedyes" "custonormal" colnames(base_teste_final)[1] "age" "sexmale" "bmi" "children"
[5] "smokeryes" "regionnorthwest" "regionsoutheast" "regionsouthwest"
[9] "marriedyes" 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)
pred_sem_regiao_log = predict(modelo_sem_regiao_log,newdata = base_teste_final)Medidas de Qualidade
Primeiro vamos guardar o valor real da variável alvo na base de teste.
y = ifelse(base_teste$custo == "normal",1,0)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 42 8
1 6 200
Accuracy : 0.9453
95% CI : (0.9099, 0.9698)
No Information Rate : 0.8125
P-Value [Acc > NIR] : 5.475e-10
Kappa : 0.8233
Mcnemar's Test P-Value : 0.7893
Sensitivity : 0.8750
Specificity : 0.9615
Pos Pred Value : 0.8400
Neg Pred Value : 0.9709
Prevalence : 0.1875
Detection Rate : 0.1641
Detection Prevalence : 0.1953
Balanced Accuracy : 0.9183
'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 42 7
1 6 201
Accuracy : 0.9492
95% CI : (0.9147, 0.9727)
No Information Rate : 0.8125
P-Value [Acc > NIR] : 1.335e-10
Kappa : 0.8347
Mcnemar's Test P-Value : 1
Sensitivity : 0.8750
Specificity : 0.9663
Pos Pred Value : 0.8571
Neg Pred Value : 0.9710
Prevalence : 0.1875
Detection Rate : 0.1641
Detection Prevalence : 0.1914
Balanced Accuracy : 0.9207
'Positive' Class : 0