Base de Seguros

um exemplo de RNA para regressão

Author

Jessica Kubrusly

O Problema

Considere a base de dados insurance.csv, composta pelas variáveis: age, sex, bmi, children, smoker, region, married, charges.

Vamos ajustar diferentes modelos de regressão para prever o custo (charge) de um cliente 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.

library(caret)
library(tidyverse)
library(neuralnet)

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…

A variável country será retirada pois para todas as observações da base temos country="usa".

table(base$country)

 usa 
1337 
base = base |> select(-country)

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 = nrow(base)
indices_treino = createDataPartition(1:N,p=0.7)[[1]]
base_treino = base[indices_treino,]
base_teste  = base[-indices_treino,]

É sempre importante salvarmos os objetos gerados para que em uma análise futura podemos começar nesta etapa e não precisaremos novamente dividir a base em treino e teste. O salvamento pode ser feito no formato RDS, que guarda os objetodos do R para futuros carregamentos.

saveRDS(base_treino,file = "base_treino_bruta.RDS")
saveRDS(base_teste,file = "base_teste_bruta.RDS")

Tratamento das NAs

Sobre o tratamento das NAs, temos a alternativa de eliminar as linhas, de imputar valores ou eliminar colunas.

Antes de qualquer análise, precisamos eliminar os indivíduos com NA na variável alvo, no caso desse problema, nas variável charge. Isso porque uma NA nessa variável não temos outra alternativa.

indices = which(is.na(base_treino$charges))
if(length(indices)>0){
  base_treino = base_treino[-indices,]  
}

O pacote naniar será útil para a análise das NAs.

library(naniar)

Nesta pacote usaremos as funções miss_var_summary e gg_miss_var para entender melhor a situação de dados faltantes na base de treino.

miss_var_summary(base_treino)
# A tibble: 8 × 3
  variable n_miss pct_miss
  <chr>     <int>    <num>
1 married     381   40.6  
2 age           1    0.107
3 bmi           1    0.107
4 sex           0    0    
5 children      0    0    
6 smoker        0    0    
7 region        0    0    
8 charges       0    0    
gg_miss_var(base_treino)

O resultado nos mostra uma grande proporção de dados faltantes na variável married, o que justifica a sua retirada.

base_treino = base_treino |> select(-married)

Repetimos o processo para ver como ficou a situação de dados faltantes na base.

miss_var_summary(base_treino)
# A tibble: 7 × 3
  variable n_miss pct_miss
  <chr>     <int>    <num>
1 age           1    0.107
2 bmi           1    0.107
3 sex           0    0    
4 children      0    0    
5 smoker        0    0    
6 region        0    0    
7 charges       0    0    
gg_miss_var(base_treino)

Com apenas 1 observação faltante em 2 variáveis, vamos eliminar as linhas que contém algum NA de forma a finalizar esta etapa.

base_treino = na.omit(base_treino)

Podemos nesse momento salvar a base de treino sem dados faltantes.

saveRDS(base_treino,file = "base_treino_sem_NA.RDS")

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.25321
attr(age_,"scaled:scale")
[1] 14.19539
bmi_ = scale(base_treino$bmi)

children_ = scale(base_treino$children)

Padronização das variáveis dependente numérica

A padronização da variável dependente tem que ser feita de forma que so valores padronizados estejam contidos na imagem da função de ativação escolhida para o problema.

No caso de usarmos a função de ativação logistica, vamos padronizar a variável resposta de forma que seus valores fiquem dentro do intervalo \([0,1]\). A operação realizada será: \[ \tilde{Y} = \dfrac{Y - \min\{Y\}}{\max\{Y\}- \min\{Y\}} \]

charges_log = (base_treino$charges - min(base_treino$charges))/(
  max(base_treino$charges)  - min(base_treino$charges)
)

Já para o caso de usarmos a função de ativação tangente hiperbólica, vamos padronizar a variável resposta de forma que seus valores fiquem dentro do intervalo \([-1,1]\). A operação realizada será: \[ \tilde{Y} = 2 \dfrac{Y - \min\{Y\}}{\max\{Y\}- \min\{Y\}} - 1 \]

charges_tan = 2*((base_treino$charges - min(base_treino$charges))/(
  max(base_treino$charges)  - min(base_treino$charges)
)) - 1

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.

Vamos fazer duas bases, uma que será usada no ajuste do modelo com função de ativação logística e outro com função de ativação tangente hiperbólica. Primeiro a base para a logística, onde a variável charge foi padronizada de forma que seus valores estão dentro de \([0,1]\).

base_treino_log = tibble(
  age = as.numeric(age_),
  bmi = as.numeric(bmi_),
  children = as.numeric(children_),
  charges = as.numeric(charges_log),
  sex = base_treino$sex,
  smoker = base_treino$smoker,
  region = base_treino$region
)
  
glimpse(base_treino_log)
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…
$ charges  <dbl> 2.514957e-01, 9.483651e-03, 5.296955e-02, 4.366851e-02, 4.190…
$ 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…
base_treino_final_log = 
model.matrix( ~ age + sex + bmi + 
                children + smoker + 
                region + charges , 
              data = base_treino_log)

base_treino_final_log = base_treino_final_log[,-1]                                
                                     
head(base_treino_final_log)
         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     charges
1               0               1 0.251495668
2               1               0 0.009483651
3               1               0 0.052969549
4               0               0 0.043668514
5               1               0 0.041908688
6               0               0 0.098181751

Agora a base para a tangente hiporbólica, onde a variável charge foi padronizada de forma que seus valores estão dentro de \([-1,1]\).

base_treino_tan = tibble(
  age = as.numeric(age_),
  bmi = as.numeric(bmi_),
  children = as.numeric(children_),
  charges = as.numeric(charges_tan),
  sex = base_treino$sex,
  smoker = base_treino$smoker,
  region = base_treino$region
)
  
glimpse(base_treino_tan)
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…
$ charges  <dbl> -0.4970087, -0.9810327, -0.8940609, -0.9126630, -0.9161826, -…
$ 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…
base_treino_final_tan = 
model.matrix( ~ age + sex + bmi + 
                children + smoker + 
                region + charges , 
              data = base_treino_tan)

base_treino_final_tan = base_treino_final_tan[,-1]                                
                                     
head(base_treino_final_tan)
         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    charges
1               0               1 -0.4970087
2               1               0 -0.9810327
3               1               0 -0.8940609
4               0               0 -0.9126630
5               1               0 -0.9161826
6               0               0 -0.8036365

Podemos nesse momento salvar a base de treino com as variáveis quantitativas padronizadas e qualitativas transformadas em indicadoras. Isso será feito considerando tanto a transformação de \(Y\) para a função de ativação logística quanto para a função de ativação tangente hiperbólica.

saveRDS(base_treino_final_log,file = "base_treino_final_log.RDS")
saveRDS(base_treino_final_tan,file = "base_treino_final_tan.RDS")

Ajuste 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 = charges ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast +  regionsouthwest, 
  data = base_treino_final_log, 
  hidden = 0,
  linear.output = FALSE)

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.043981222  0.359741255  0.005402707  0.332734078  0.048987691
[6]  2.114865732 -0.045748489 -0.103676165 -0.095012211

Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:

head(modelo_completo_log$net.result[[1]])
        [,1]
1 0.32438099
2 0.07442389
3 0.09740664
4 0.08227724
5 0.06481513
6 0.09722668

Faç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])
[1] -2.043981222  0.359741255  0.005402707  0.332734078  0.048987691
[6]  2.114865732 -0.045748489 -0.103676165 -0.095012211
(x = as.numeric(modelo_completo_log$data[1,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest")]))
[1] -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000  0.0000000  0.0000000
[8]  1.0000000
(x = c(1,x))
[1]  1.0000000 -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000  0.0000000
[8]  0.0000000  1.0000000
soma = sum(x*pesos)
1/(1+exp(-soma))
[1] 0.324381
modelo_completo_log$net.result[[1]][1]
[1] 0.324381

Modelo Completo função de ativação tengente hiperbólica

modelo_completo_tan = neuralnet(
  formula = charges ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast +  regionsouthwest, 
  data = base_treino_final_tan, 
  hidden = 0,
  linear.output = FALSE,
  act.fct = "tanh")

Podemos visualizar o modelo criado com a função plot:

plot(modelo_completo_tan,rep = "best")

Podemos acessar os valores estimados pelos pesos sinápticos da seguinte forma:

modelo_completo_tan$weights[[1]][[1]][,1]
[1] -1.023708313  0.179971658  0.003456346  0.166348883  0.024421306
[6]  1.057889456 -0.021748421 -0.050866209 -0.046428991

Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:

head(modelo_completo_tan$net.result[[1]])
        [,1]
1 -0.3514598
2 -0.8511922
3 -0.8052547
4 -0.8353900
5 -0.8705449
6 -0.8057975

Faça na mão a estimativa para a primeira observação da base e verifique se está correta:

(pesos = modelo_completo_tan$weights[[1]][[1]][,1])
[1] -1.023708313  0.179971658  0.003456346  0.166348883  0.024421306
[6]  1.057889456 -0.021748421 -0.050866209 -0.046428991
(x = as.numeric(modelo_completo_tan$data[1,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest")]))
[1] -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000  0.0000000  0.0000000
[8]  1.0000000
(x = c(1,x))
[1]  1.0000000 -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000  0.0000000
[8]  0.0000000  1.0000000
soma = sum(x*pesos)
(exp(2*soma)-1)/(exp(2*soma)+1)
[1] -0.3514598
modelo_completo_tan$net.result[[1]][1]
[1] -0.3514598

Modelo sem Região com função de ativação logística

modelo_sem_regiao_log = neuralnet(
  formula = charges ~ age + sexmale + bmi + children + smokeryes , 
  data = base_treino_final_log, 
  hidden = 0,
  linear.output = FALSE)

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] -2.105064188  0.362891866  0.001986733  0.322036339  0.047292124
[6]  2.112519398

Podemos acessar os valores estimados para os dados de treinamento da seguinte forma:

head(modelo_sem_regiao_log$net.result[[1]])
        [,1]
1 0.33176798
2 0.07646891
3 0.10007832
4 0.08110472
5 0.06797063
6 0.09605654

Faç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])
[1] -2.105064188  0.362891866  0.001986733  0.322036339  0.047292124
[6]  2.112519398
(x = as.numeric(modelo_sem_regiao_log$data[1,c("age","sexmale","bmi","children","smokeryes")]))
[1] -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000
(x = c(1,x))
[1]  1.0000000 -1.4267453  0.0000000 -0.4542083 -0.9225408  1.0000000
soma = sum(x*pesos)
(1)/(1+exp(-soma))
[1] 0.331768
modelo_sem_regiao_log$net.result[[1]][1]
[1] 0.331768

Medidas de Qualidade na base de treino - problemas de Regressão

Vamos comparar o desempenho dos três modelos de RNA ajustados para prever o custo 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_completo_tan_ = modelo_completo_tan$net.result[[1]]
y_sem_regiao_log_ = modelo_sem_regiao_log$net.result[[1]]

Repare que o valor previsto pelo modelo é para a variável resposta transformada. Mas para calcular os erros de previsão precisamos da previsão para a variável na unidade original e para isso precisaremos “desfazer” a transformação.

Para desfazer a transformação feita para a função de ativação logística:

\[ Y = \tilde{Y} \times (\max\{Y\}- \min\{Y\}) + \min\{Y\} \]

por isso precisamos ter guadrado os valores usados para transformação.

min = min(base_treino$charges)
max = max(base_treino$charges)
y_completo_log = y_completo_log_ * (max - min) + min
y_sem_regiao_log = y_sem_regiao_log_ * (max - min) + min

Para desfazer a transformação feita para a função de ativação tangente hiperbólica:

\[ Y = \dfrac{\tilde{Y} + 1}{2} \times (\max\{Y\}- \min\{Y\}) + \min\{Y\} \]

y_completo_tan = ((y_completo_tan_ + 1)/2) * (max - min) + min

Para facilitar as contas que seguem, vamos primeiro definir o erro cometido em cada modelo para cada indivíduo.

erro_completo_log = y_completo_log - base_treino$charges
erro_completo_tan = y_completo_tan - base_treino$charges
erro_sem_regiao_log = y_sem_regiao_log - base_treino$charges

O gráfico dos erros pode nos indicar se alguma observaçao da base é está muito fora do padrão e deve ser retirada da base de treinamento.

plot(erro_completo_log)

plot(erro_completo_tan)

plot(erro_sem_regiao_log)

Soma dos Erros Quadráticos

A Soma dos Erros Quadráticos (SSE - Sum of Squared Errors) é definida por

\[ SSE = \sum_{i=1}^N (y_i - \hat{y}_i)^2 \]

No R podemos calculá-la assim:

(sse_completo_log = sum((erro_completo_log)^2))
[1] 31926894351
(sse_completo_tan = sum((erro_completo_tan)^2))
[1] 31926709505
(sse_sem_regiao_log = sum((erro_sem_regiao_log)^2))
[1] 32044644087

Erro Quadrático Médio

O Erro Quadrático Médio (Mean Squared Error - MSE) é definido por

\[ MSE = \dfrac{1}{N} \sum_{i=1}^N (y_i - \hat{y}_i)^2 \]

No R podemos calculá-la assim:

n = nrow(base_treino)
(mse_completo_log = sse_completo_log/n)
[1] 34109930
(mse_completo_tan = sse_completo_tan/n)
[1] 34109732
(mse_sem_regiao_log = sse_sem_regiao_log/n)
[1] 34235731

Coeficiente de Determinação \(R^2\)

O Coeficiente de Determinação \(R^2\) é definido por:

\[ R^2 = 1 - \dfrac{SSE}{SST} = 1 - \dfrac{\sum_{i=1}^N (y_i - \hat{y}_i)^2 }{\sum_{i=1}^N (y_i - \bar{y})^2 } \quad \text{ sendo, }\quad SST = \sum_{i=1}^N (y_i - \bar{y})^2 \]

sst = sum((mean(base_treino$charges) - base_treino$charges)^2)
(R2_completo_log = 1 - sse_completo_log/sst)
[1] 0.7614635
(R2_completo_tan = 1 - sse_completo_tan/sst)
[1] 0.7614648
(R2_sem_regiao_log = 1 - sse_sem_regiao_log/sst)
[1] 0.7605837

Medidas de Qualidade na base de teste - problemas de Regressão

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.

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…
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_children

Agora vamos criar a base de teste transformada (sem a variável resposta - charge) e 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_tan)
[1] "age"             "sexmale"         "bmi"             "children"       
[5] "smokeryes"       "regionnorthwest" "regionsoutheast" "regionsouthwest"
[9] "charges"        
colnames(base_treino_final_log)
[1] "age"             "sexmale"         "bmi"             "children"       
[5] "smokeryes"       "regionnorthwest" "regionsoutheast" "regionsouthwest"
[9] "charges"        
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_completo_tan_ = predict(modelo_completo_tan,newdata = base_teste_final)[,1]
pred_sem_regiao_log_ = predict(modelo_sem_regiao_log,newdata = base_teste_final)[,1]

E podemos retornar a previsão da variável resposta para a unidade original.

pred_completo_log = pred_completo_log_ * (max - min) + min
pred_completo_tan = ((pred_completo_tan_ + 1)/2)* (max - min) + min
pred_sem_regiao_log = pred_sem_regiao_log_ * (max - min) + min

Medidas de Qualidade

O processo será o mesmo realizado na base de treino.

Primeiro vamos definir o erro cometido em cada modelo para cada indivíduo.

erro_completo_log_teste = pred_completo_log - base_teste$charges
erro_completo_tan_teste = pred_completo_tan - base_teste$charges
erro_sem_regiao_log_teste = pred_sem_regiao_log -  base_teste$charges
plot(erro_completo_log_teste)

plot(erro_completo_tan_teste)

plot(erro_sem_regiao_log_teste)

(sse_completo_log_teste = sum((erro_completo_log_teste)^2))
[1] 14101692421
(sse_completo_tan_teste = sum((erro_completo_tan_teste)^2))
[1] 14098665113
(sse_sem_regiao_log_teste = sum((erro_sem_regiao_log_teste)^2))
[1] 14200545164
n = nrow(base_teste)
(mse_completo_log_teste = sse_completo_log_teste/n)
[1] 35520636
(mse_completo_tan_teste = sse_completo_tan_teste/n)
[1] 35513010
(mse_sem_regiao_log_teste = sse_sem_regiao_log_teste/n)
[1] 35769635

Atenção no código a seguir, a média usada no cálculo do SST será sempre a da base de treino.

sst_teste = sum((mean(base_treino$charges) - base_teste$charges)^2)
(R2_completo_log_teste = 1 - sse_completo_log_teste/sst_teste)
[1] 0.7631603
(R2_completo_tan_teste = 1 - sse_completo_tan_teste/sst_teste)
[1] 0.7632111
(R2_sem_regiao_log_teste = 1 - sse_sem_regiao_log_teste/sst_teste)
[1] 0.7615

Comparação entre treino e teste:

valores = 
  matrix(
    c(R2_completo_log,R2_completo_tan, R2_sem_regiao_log,
      R2_completo_log_teste,R2_completo_tan_teste,R2_sem_regiao_log_teste), nrow=2, byrow = T)
valores
          [,1]      [,2]      [,3]
[1,] 0.7614635 0.7614648 0.7605837
[2,] 0.7631603 0.7632111 0.7615000
barplot(valores,beside = T, ylim = c(0,1),col=c("tomato","blue3"),
  legend.text = c("A","B"), names.arg = c("completo log","completo tan","sem regiao"))