Explorando um navio afundado
Autores: Filipe Coelho e Leonardo Ferreira
20 de Junho de 2019
1 INTRODUÇÃO
Iremos apresentar nosso trabalho utilizando a base de dados do Titanic. É uma base famosa e de fácil compreensão. Utilizaremos três técnicas de previsão que foram apresentados durante o curso de Data Science: * ÁRVORE DE DECISÃO * PROBIT * LOGIT
Também faremos análise gráficas para ilustrar e mensurar o trabalho feito.
O objetivo do trabalho é utilizar as variáveis de características dos passageiros para prever as mortes do acidente.
Estruturamos o script em basicamente quatro etapas:
- Carregamento de pacotes e Importação de Dados
- ETL e input de Missing values
- Previsão
- Conclusão
2 PACOTES E IMPORTAÇÃO DA BASE DE DADOS
Comando que verifica se o pacote está instalado. Se não estiver, instala e carrega.
if(!require(rpart)) {install.packages("rpart")}
require(rpart)
if(!require(rattle)) {install.packages("rattle")}
require(rattle)
if(!require(rpart.plot)) {install.packages("rpart.plot")}
require(rpart.plot)
if(!require(readxl)) {install.packages("readxl")}
require(readxl)
if(!require(readr)) {install.packages("readr")}
require(readr)
if(!require(caret)) {install.packages("caret")}
require(caret)
if(!require(tidyverse)) {install.packages("tidyverse")}
require(tidyverse)
if(!require(ggplot2)) {install.packages("ggplot2")}
require(ggplot2)
if(!require(ggthemes)) {install.packages("ggthemes")}
require(ggthemes)
if(!require(mice)) {install.packages("mice")}
require(mice)Importando Base de Dados
titanic <- read_excel("titanic.xls", col_types = c("numeric",
"numeric", "text", "text", "numeric",
"numeric", "numeric", "text", "numeric",
"text", "text", "text", "numeric", "text",
"numeric"))Agora que temos nossas varáveis, sabemos que temos 1309 observações e 12 variáveis Abaixo segue uma breve descrição de cada uma:
| Nome | Descrição |
|---|---|
| Survived | Sobreviveu (1) ou morreu (0) |
| Pclass | Classe do passageiro |
| Name | Nome do passageiro |
| Sex | Sexo do passageiro |
| Age | Idade do passageiro |
| SibSp | Número de irmãos/esposo(a) a bordo |
| Parch | Número de pais/filhos(as) a bordo |
| Ticket | Número do ticket |
| Fare | Preço do ticket |
| Cabin | Cabine |
| Home.dest | Destino do passageiro |
| Embarked | Porto de embarque |
Vamos analisar as estatísticas descritivas e características das variáveis conforme elas estavam originalmente no banco de dados (sem nenhum tratamento):
summary(titanic)## pclass survived name sex
## Min. :1.000 Min. :0.000 Length:1309 Length:1309
## 1st Qu.:2.000 1st Qu.:0.000 Class :character Class :character
## Median :3.000 Median :0.000 Mode :character Mode :character
## Mean :2.295 Mean :0.382
## 3rd Qu.:3.000 3rd Qu.:1.000
## Max. :3.000 Max. :1.000
##
## age sibsp parch ticket
## Min. : 0.17 Min. :0.0000 Min. :0.000 Length:1309
## 1st Qu.:21.00 1st Qu.:0.0000 1st Qu.:0.000 Class :character
## Median :28.00 Median :0.0000 Median :0.000 Mode :character
## Mean :29.56 Mean :0.4992 Mean :0.385
## 3rd Qu.:38.00 3rd Qu.:1.0000 3rd Qu.:0.000
## Max. :80.00 Max. :8.0000 Max. :9.000
## NA's :9 NA's :1
## fare cabin embarked
## Min. : 0 Length:1309 Length:1309
## 1st Qu.: 14 Class :character Class :character
## Median : 74 Mode :character Mode :character
## Mean : 130656
## 3rd Qu.: 78958
## Max. :5123292
## NA's :1
## boat body home.dest PassengerId
## Length:1309 Min. : 1.0 Length:1309 Min. : 1.0
## Class :character 1st Qu.: 72.0 Class :character 1st Qu.: 327.0
## Mode :character Median :155.0 Mode :character Median : 654.0
## Mean :160.8 Mean : 654.4
## 3rd Qu.:256.0 3rd Qu.: 982.0
## Max. :328.0 Max. :1309.0
## NA's :1188str(titanic)## Classes 'tbl_df', 'tbl' and 'data.frame': 1309 obs. of 15 variables:
## $ pclass : num 3 1 3 1 3 3 1 3 3 2 ...
## $ survived : num 0 1 1 1 0 0 0 0 1 1 ...
## $ name : chr "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
## $ sex : chr "male" "female" "female" "female" ...
## $ age : num 22 38 16 35 35 22 54 2 27 14 ...
## $ sibsp : num 1 1 0 1 0 0 0 3 0 1 ...
## $ parch : num 0 0 0 0 0 0 0 1 2 0 ...
## $ ticket : chr "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...
## $ fare : num 7.25 7.13e+05 7.92e+03 5.31e+01 8.05 ...
## $ cabin : chr NA "C85" NA "C123" ...
## $ embarked : chr "S" "C" "S" "S" ...
## $ boat : chr NA "4" NA "D" ...
## $ body : num NA NA NA NA NA NA 175 4 NA NA ...
## $ home.dest : chr "UK" "USA" "FINLAND" "USA" ...
## $ PassengerId: num 1 2 3 4 5 6 7 8 9 10 ...Verifica-se que várias variáveis categóricas não estão classificadas como factor. Além disso, algumas delas apresentam valores em branco (NA). Dessa forma, será realizado abaixo o tratamento necessário dessas variáveis.
3 ETL e INPUT DE MISSING VALUE
3.1 Tratamento de dados
Transformamos as variáveis categóticas em factor, substituimos os NA’s na coluna Fare pela mediana das células com valor, tratamos a variável name para extrair dados que sejam úteis na regressão ou na geração da nova variável (o título do passageiro: Title) e excluímos algumas variveis que não contribuiriam para a previsão, como o número do ticket.
combi <- titanic
combi$name <- as.character(combi$name)
strsplit(combi$name[1], split='[,.]')## [[1]]
## [1] "Braund" " Mr" " Owen Harris"strsplit(combi$name[1], split='[,.]')[[1]]## [1] "Braund" " Mr" " Owen Harris"strsplit(combi$name[1], split='[,.]')[[1]][2]## [1] " Mr"combi$Title <- sapply(combi$name, FUN=function(x) {strsplit(x, split='[,.]')[[1]][2]})
combi$Title <- sub(' ', '', combi$Title)
combi$Title[combi$PassengerId == 797] <- 'Mrs' # ajustdo para Mrs.porque era mulher
combi$Title[combi$Title %in% c('Lady', 'the Countess', 'Mlle', 'Mee', 'Ms', 'Mme')] <- 'Miss'
combi$Title[combi$Title %in% c('Capt', 'Don', 'Major', 'Sir', 'Col', 'Jonkheer', 'Rev', 'Dr')] <- 'Mr'
combi$Title[combi$Title %in% c('Dona')] <- 'Mrs'
combi$Title <- as.factor(combi$Title)
unique(combi$Title)## [1] Mr Mrs Miss Master
## Levels: Master Miss Mr Mrscombi$embarked <- factor(combi$embarked)
combi$pclass <- as.factor(combi$pclass)
combi$survived <- as.factor(combi$survived)
combi$sex <- as.factor(combi$sex)
combi$fare[is.na(combi$fare)] <- median(combi$fare, na.rm=TRUE)
excluir <- c("ticket", "cabin", "boat", "body", "home.dest")
combi <- combi[,!(names(combi)%in% excluir)]Ser rico compensa na hora da sobrevivência? Vamos plotar dois gráficos para podermos entendermos melhor a relação das variáveis Pclass e Title com a sobrevivência.
# Relação entre title e survived
ggplot(combi, aes(x = Title, fill = factor(survived))) +
geom_bar(stat='count', position='dodge') +
labs(x = 'Title') 
O gráfico acima apresenta a quantidade de pessoas que sobreviveram e que morreram em cada um dos grupos da variável Title. Verificamos que mais mulheres sobreviveram que morreram, tanto entre as solteiras quanto entre as casadas. Entre os homens, a grande maioria morreu. Esse gráfico já antecipa que o gênero será uma variável importante na previsão
# Relação entre pclass e survived
ggplot(combi, aes(x = pclass, fill = factor(survived))) +
geom_bar(stat='count', position='dodge') +
labs(x = 'Classe')
No gráfico acima, analisamos a quantidade de pessoas que osbreviveram em cada uma das classes de passageiros. Conclui-se que dentre os passageiros da primeira classe, houve mais sobreviventes que mortos. Na segunda classe, houve mais mortos que sobrevientes. Na terceira classe, aproximadamente 3/4 dos passageiros morreram. Dessa forma, é possível antecipar que a classe do passageiro será uma variável improtante na previsão.
Famílias afundam ou nadam juntas?
Depois de dividir os títulos de cada passageiros, vamos avançar e fazer variáveis familiares. Primeiro, vamos criar a variável family_size baseada no número de irmãos/conjuges e no número de filhos/parentes a bordo.
combi$Surname <- sapply(combi$name,
function(x) strsplit(x, split = '[,.]')[[1]][1])
combi$family_size <- combi$sibsp + combi$parch + 1
# Substituindo os NA's na colune family Size pela mediana das células com valor.
combi$family_size[is.na(combi$family_size)] <- median(combi$family_size, na.rm=TRUE)Vamos plotar esta nova variável family_size para entender melhor como ela pode estar relacionada à sobrevivência.
ggplot(combi[1:891,], aes(x = family_size, fill = factor(survived))) +
geom_bar(stat='count', position='dodge') +
scale_x_continuous(breaks=c(1:11)) +
labs(x = 'Family Size') +
theme_few()
Inferimos do gráfico no grupo de solteiros houve mais mortos que sobreviventes. O mesmo aconteceu com famílias com 5 ou mais pessoas. Enquanto isso, nas famílias com 2 a 4 pessoas, houve mais sobreviventes que mortos. Com isso, decidimos criar três novas variáveis categóricas: Solteiros, pequeno (de 1 à 4 familiares) e grande (5 ou mais familiares a bordo) para melhor visualização. O nome da variável será: family_size_discreto
# Tamanho da família discreto
combi$family_size_discreto[combi$family_size == 1] <- 'solteiro'
combi$family_size_discreto[combi$family_size < 5 & combi$family_size > 1] <- 'pequeno'
combi$family_size_discreto[combi$family_size > 4] <- 'grande'
combi$family_size_discreto <- factor(combi$family_size_discreto)
# Mosaico com tamanho familiar por sobrevivente
mosaicplot(table(combi$family_size, combi$survived), main='Tamanho da Família por nº de sobreviventes', shade=TRUE)
O mosaico acima ajuda a reforçar a ideia de que famílias pequenas possuem mais chance e sobreviver (porém não os solteiros).
3.2 Prevendo idades com NA
Fazemos a previsão das idades vazias de cada passageiro pelo método de árvore de decisão, utilizando as outras variáveis da amostra. Foi utilizado o method=“anova” que é um método de regression tree, pois estamos prevendo uma variável contínua.
predicted_age <- rpart(age ~ pclass + sex + sibsp + parch + fare + embarked + Title + family_size,
data=combi[!is.na(combi$age),], method="anova")
combi$age[is.na(combi$age)] <- predict(predicted_age, combi[is.na(combi$age),])Vamos analisar o gráfico que apresenta a distribuição das idades dos grupos de sobreviventes e mortos:
#relação entre age e survived
ggplot(combi, aes(x = age)) +
geom_histogram(aes(fill = survived), position = "dodge", stat = "density") +
ylab("Density") +
xlab("Age") +
ggtitle("Distribuição de idade")
É possível ver nesse gráfico, que entre as crianças, houve mais sobrevivente que mortos.
O gráfico seguinte mostra a quantidade de sobreviventes e mortos por idade, porém separados por sexo:
# Relação idade x sobrevivência
ggplot(combi, aes(age, fill = factor(survived))) +
geom_histogram() +
facet_grid(.~sex) +
theme_few()
Além de reforçar que as crianças possuem boas chances de sobreviver, também reforça que os homens tem mais chances de morrer e que o grupo mais numeroso a bordo era o de homens entre 20 e 50 anos, tendo sobrevivido a menor parte deste grupo.
Agora que temos as idades completas, criaremos mais duas colunas de variáveis categóricas: child e mother. Se o passageiro tiver menos de 18 anos, será considerado crianca, caso contrário adulto. Se o passageiro tiver mais de 18 anos, for do sexo feminino, tiver mais de um membro familiar a bordo e tiver seu Title diferente de Miss (senhorita), será considerado uma mãe.
# Adicionando a variável **Child**
combi$child[combi$age < 18] <- 'crianca'
combi$child[combi$age >= 18] <- 'adulto'
# Tabela com quantidades:
table(combi$child, combi$survived)##
## 0 1
## adulto 716 415
## crianca 93 85Parece que ser criança não é tão ruim assim, mas não é garantia de sobreviVência também. Como sabemos (pelo filme), mulheres e crianças foram colocados primeiros nos botes salva vidas. Talvez as mães tenham se salvado também, veremos a seguir:
# Adicionando a variável **Mother**
combi$Mother <- 'nao_mae'
combi$Mother[combi$sex == 'female' & combi$parch > 0 & combi$age > 18 & combi$Title != 'Miss'] <- 'mae'
combi$child <- factor(combi$child)
combi$Mother <- factor(combi$Mother)
# Tabela com quantidades:
table(combi$Mother, combi$survived)##
## 0 1
## mae 21 64
## nao_mae 788 436De fato, ser mãe também ajudou a aumentar as chances de sobreviver.
Depois de todos os ajustes, vamos analisar novamente as estatísticas descritivas e características das variáveis:
summary(combi)## pclass survived name sex age
## 1:323 0:809 Length:1309 female:466 Min. : 0.17
## 2:277 1:500 Class :character male :843 1st Qu.:21.00
## 3:709 Mode :character Median :28.00
## Mean :29.61
## 3rd Qu.:38.00
## Max. :80.00
##
## sibsp parch fare embarked
## Min. :0.0000 Min. :0.000 Min. : 0 C:270
## 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.: 14 Q:123
## Median :0.0000 Median :0.000 Median : 74 S:916
## Mean :0.4992 Mean :0.385 Mean : 130556
## 3rd Qu.:1.0000 3rd Qu.:0.000 3rd Qu.: 78958
## Max. :8.0000 Max. :9.000 Max. :5123292
## NA's :1
## PassengerId Title Surname family_size
## Min. : 1.0 Master: 61 Length:1309 Min. : 1.000
## 1st Qu.: 327.0 Miss :267 Class :character 1st Qu.: 1.000
## Median : 654.0 Mr :782 Mode :character Median : 1.000
## Mean : 654.4 Mrs :199 Mean : 1.884
## 3rd Qu.: 982.0 3rd Qu.: 2.000
## Max. :1309.0 Max. :11.000
##
## family_size_discreto child Mother
## grande : 82 adulto :1131 mae : 85
## pequeno :437 crianca: 178 nao_mae:1224
## solteiro:790
##
##
##
## str(combi)## Classes 'tbl_df', 'tbl' and 'data.frame': 1309 obs. of 16 variables:
## $ pclass : Factor w/ 3 levels "1","2","3": 3 1 3 1 3 3 1 3 3 2 ...
## $ survived : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 1 2 2 ...
## $ name : chr "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
## $ sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
## $ age : num 22 38 16 35 35 22 54 2 27 14 ...
## $ sibsp : num 1 1 0 1 0 0 0 3 0 1 ...
## $ parch : num 0 0 0 0 0 0 0 1 2 0 ...
## $ fare : num 7.25 7.13e+05 7.92e+03 5.31e+01 8.05 ...
## $ embarked : Factor w/ 3 levels "C","Q","S": 3 1 3 3 3 2 3 3 3 1 ...
## $ PassengerId : num 1 2 3 4 5 6 7 8 9 10 ...
## $ Title : Factor w/ 4 levels "Master","Miss",..: 3 4 2 4 3 3 3 1 4 4 ...
## $ Surname : chr "Braund" "Cumings" "Heikkinen" "Futrelle" ...
## $ family_size : num 2 2 1 2 1 1 1 5 3 2 ...
## $ family_size_discreto: Factor w/ 3 levels "grande","pequeno",..: 2 2 3 2 3 3 3 1 2 2 ...
## $ child : Factor w/ 2 levels "adulto","crianca": 1 1 2 1 1 1 1 2 1 2 ...
## $ Mother : Factor w/ 2 levels "mae","nao_mae": 2 2 2 2 2 2 2 2 1 2 ...Alguns destaques da estatística descritiva desses dados ajustados: - A maioria das pessoas era da terceira classe; - A maioria morreu; - Os homens eram quase o dobro das mulheres; - A mediana da idade era 28 anos; - A maioria viajava sozinho; - Havia poucas crianças e poucas mães.
4 CRIANDO MODELO DE PREVISÃO
4.1 Preparando base de dados para Cross Validation
No K-Fold Cross Validation, você divide normalmente entre 5–10 grupos o seu conjunto de dados. 1 grupo é o grupo de teste e k-1 grupos juntos são o grupo de treino. O modelo “roda” k vezes, alternando grupo de teste, de formaque seja possível testar em todo o conjunto de dados. Em cada rodada, é gerada uma medida de erro. Neste trabalho, optamos usar Matriz de Confusão. A partir dela, extraímos uma medida de acurácia, que é o percentual de acerto da prevsião. A medida final do modelo é a média desta acurácia em cada uma das k rodadas. Utilizaremos essa técnica em todos os modelos.
k_folds <- 5
folds <- createFolds(combi$survived, k = k_folds, list = FALSE)
useful_cols = setdiff(colnames(combi), c("name", "PassengerId", "Surname", "family_size"))
accuracy <- numeric(k_folds)4.2 Criando o modelo de Árvore de Decisão
Primeiramente criamos uma àrvore de Decisão com o CP padrão (0.01)
accuracy_tree <- numeric(k_folds)
for (i in 1:k_folds) {
data_train_tree = combi[folds != i, useful_cols]
data_test_tree = combi[folds == i, useful_cols]
Tree <- rpart(survived ~
pclass+ sex+ age+
fare+
embarked+
Title+
family_size_discreto+
child+
Mother, data = data_train_tree)
cast_tree <- predict(Tree, newdata = data_test_tree)
prediction_tree <- table(data_test_tree$survived, cast_tree[ ,"1"] > 0.5)
TP_tree <- prediction_tree["1", "TRUE"]
FP_tree <- prediction_tree["0", "TRUE"]
TN_tree <- prediction_tree["0", "FALSE"]
FN_tree <- prediction_tree["1", "FALSE"]
accuracy_tree[i] = (TP_tree + TN_tree)/(TP_tree + TN_tree + FP_tree + FN_tree)}
paste("Acurácia média", toString(k_folds), "folds", toString(round(mean(accuracy_tree), 4)))## [1] "Acurácia média 5 folds 0.8113"Visualizando a árvore de decisão
fancyRpartPlot(Tree)
Vamos agora buscar o melhor CP (como menor erro de validação):
#Escolha valores razoáveis
possible_values <- seq(from = 0.0001, to = 0.1, length.out = 60)
cp_error <- numeric(length(possible_values))
#Escolha a quantidade de dados que serão alocado para teste
p <- 0.2
test_idx <- sample(1:nrow(combi), size = floor(p*nrow(combi)))
data_test_cp <- combi[test_idx, useful_cols]
data_train_val_cp <- combi[-test_idx, useful_cols]
iteration <- 1
for (val in possible_values) {
validation_errors <- numeric(k_folds)
# Para cada valor do parametro cp teremos a seguinte rotina
for (run in 1:k_folds) {
# A cada rodada de validação, separamos em dados de treino e de validação
df_train <- data_train_val_cp[folds != run, useful_cols]
df_val <- data_train_val_cp[folds == run, useful_cols]
# Treinamos o modelo
tree_cp <- rpart(survived ~
pclass+ sex+ age+
fare+
embarked+
Title+
family_size_discreto+
child+
Mother, data = df_train, cp = val)
# Fazemos a previsão no conjunto de validação
pred_cp <- predict(tree_cp, newdata = df_val, type = "class")
# Computamos uma métrica de avaliação do erro. Exemplo: porcentagem de acerto
right <- pred_cp == df_val$survived
validation_errors[run] = length(right[right == FALSE])/nrow(df_val)}
# Depois de computar a métrica de erro em todos os folds, tomamos uma média e guardamos
cp_error[iteration] = mean(validation_errors)
iteration <- iteration + 1}
# Plotar o erro de validação
qplot(possible_values, cp_error, geom = "line")
O gráfico acima nos ajuda a ver a relação entre os CPs e os respectivos erros de validação.
Agora estimamos o modelo na união dos dados de treino e validação com o melhor CP:
best <- possible_values[which.min(cp_error)]
best## [1] 0.02211186best_tree <- rpart(survived ~
pclass+ sex+ age+
fare+
embarked+
Title+
family_size_discreto+
child+
Mother, data = data_train_val_cp, cp = best)
fancyRpartPlot(best_tree)
Estimamos o erro do modelo fora da amostra computando o erro no conjunto de teste
pred_best_tree <- predict(best_tree, newdata = data_test_cp, type = "class")
right_best_tree <- pred_best_tree == data_test_cp$survived
final_error <- length(right_best_tree[right_best_tree == FALSE])/nrow(data_test_cp)
paste("Erro final (%)", round(100 * final_error, 4))## [1] "Erro final (%) 18.0077"Finalmente, vamos estimar uma nova Árvore de Decisão com o melhor CP, repetindo o procedimento do k-fold cross validation:
accuracy_tree2 <- numeric(k_folds)
for (i in 1:k_folds) {
data_train_tree2 = combi[folds != i, useful_cols]
data_test_tree2 = combi[folds == i, useful_cols]
Tree2 <- rpart(survived ~
pclass+ sex+ age+
fare+
embarked+
Title+
family_size_discreto+
child+
Mother, data = data_train_tree2, cp=best)
cast_tree2 <- predict(Tree2, newdata = data_test_tree2)
prediction_tree2 <- table(data_test_tree2$survived, cast_tree2[ ,"1"] > 0.5)
TP_tree2 <- prediction_tree2["1", "TRUE"]
FP_tree2 <- prediction_tree2["0", "TRUE"]
TN_tree2 <- prediction_tree2["0", "FALSE"]
FN_tree2 <- prediction_tree2["1", "FALSE"]
accuracy_tree2[i] = (TP_tree2 + TN_tree2)/(TP_tree2 + TN_tree2 + FP_tree2 + FN_tree2)}
paste("A Acurácia média em", toString(k_folds), "folds é", toString(round(mean(accuracy_tree2), 4)))## [1] "A Acurácia média em 5 folds é 0.8143"fancyRpartPlot(Tree2) 
Essa representação gráfica da árvore de decisão é de fácil interpretação. A primeira caixinha representa 100% dos dados. Dentro dela, na primeira linha, você lê o zero (sinônimo de “morreu”) porque, no grupo com todos os passageiros, a maior probabilidade é de morrer. Isso vem explicitado na segunda linha da caixinha, que traz as probabilidades: 62% de chances de morrer e 38% de chances de sobreviver. Isso também se reflete na cor da caixinha (verde para quem tem mais chances de morrer e azul para quem tem mais chances de sobreviver). A terceira linha da caixinha traz o percentual que aquela caixinha representa do total de passageiros (nesse caso, 100%). Abaixo dela, vem a condição para abrir o nó em dois caminhos. Nesse caso, a pergunta é se o título do passageiro é igual a “Mr.”. Caso positivo, segue alinha da esquerda e, caso negativo, segue a linha da direita. A árvore deve ser lida seguindo essa lógica até o final. O resultado dessa árvore condiz com as percepções geradas a partir da leitura dos gráficos e das estatísticas descritivas
4.3 Modelo LOGIT:
Utilizando o modelo LOGIT para prever os sobreviventes:
for (i in 1:k_folds) {
data_train = combi[folds != i, useful_cols]
data_test = combi[folds == i, useful_cols]
Logit <- glm(survived ~ I(pclass==1)+I(pclass==2)+
I(sex=="female")+
age+
fare+
I(embarked=="C")+I(embarked=="Q")+
I(Title=="Master")+I(Title=="Miss")+I(Title=="Mr")+
I(family_size_discreto == "solteiro")+I(family_size_discreto == "grande")+
I(child=="crianca")+
I(Mother== "mae"),
data = data_train, family=binomial())
cast <- predict(Logit, newdata = data_test, type = "response")
prediction <- table(data_test$survived, cast > 0.5)
TP <- prediction["1", "TRUE"]
FP <- prediction["0", "TRUE"]
TN <- prediction["0", "FALSE"]
FN <- prediction["1", "FALSE"]
accuracy[i] = (TP + TN)/(TP + TN + FP + FN)}
paste("A Acurácia média em", toString(k_folds), "folds é", toString(round(mean(accuracy), 4)))## [1] "A Acurácia média em 5 folds é 0.8182"Abaixo, segue o resumo do modelo Logit:
summary(Logit)##
## Call:
## glm(formula = survived ~ I(pclass == 1) + I(pclass == 2) + I(sex ==
## "female") + age + fare + I(embarked == "C") + I(embarked ==
## "Q") + I(Title == "Master") + I(Title == "Miss") + I(Title ==
## "Mr") + I(family_size_discreto == "solteiro") + I(family_size_discreto ==
## "grande") + I(child == "crianca") + I(Mother == "mae"), family = binomial(),
## data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3056 -0.5940 -0.4043 0.5777 3.0727
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) -1.758e+00 3.428e-01 -5.129
## I(pclass == 1)TRUE 2.015e+00 2.655e-01 7.590
## I(pclass == 2)TRUE 9.624e-01 2.286e-01 4.210
## I(sex == "female")TRUE 3.282e+00 3.370e-01 9.739
## age -2.977e-02 8.910e-03 -3.341
## fare 2.398e-07 2.309e-07 1.038
## I(embarked == "C")TRUE 5.118e-01 2.279e-01 2.245
## I(embarked == "Q")TRUE 1.513e-01 3.046e-01 0.497
## I(Title == "Master")TRUE 2.488e+00 4.815e-01 5.168
## I(Title == "Miss")TRUE -5.983e-01 3.734e-01 -1.603
## I(Title == "Mr")TRUE NA NA NA
## I(family_size_discreto == "solteiro")TRUE 1.283e-01 2.109e-01 0.608
## I(family_size_discreto == "grande")TRUE -2.096e+00 4.105e-01 -5.106
## I(child == "crianca")TRUE -3.026e-01 3.540e-01 -0.855
## I(Mother == "mae")TRUE 8.254e-01 5.129e-01 1.609
## Pr(>|z|)
## (Intercept) 2.92e-07 ***
## I(pclass == 1)TRUE 3.20e-14 ***
## I(pclass == 2)TRUE 2.55e-05 ***
## I(sex == "female")TRUE < 2e-16 ***
## age 0.000834 ***
## fare 0.299130
## I(embarked == "C")TRUE 0.024747 *
## I(embarked == "Q")TRUE 0.619454
## I(Title == "Master")TRUE 2.37e-07 ***
## I(Title == "Miss")TRUE 0.109017
## I(Title == "Mr")TRUE NA
## I(family_size_discreto == "solteiro")TRUE 0.543027
## I(family_size_discreto == "grande")TRUE 3.30e-07 ***
## I(child == "crianca")TRUE 0.392593
## I(Mother == "mae")TRUE 0.107592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1392.63 on 1046 degrees of freedom
## Residual deviance: 887.08 on 1033 degrees of freedom
## AIC: 915.08
##
## Number of Fisher Scoring iterations: 54.4 Modelo PROBIT:
AGora, vamos utilizar o modelo PROBIT para prever os sobreviventes:
accuracy2 <- numeric(k_folds)
for (i in 1:k_folds) {
data_train2 = combi[folds != i, useful_cols]
data_test2 = combi[folds == i, useful_cols]
Probit <- glm(survived ~ I(pclass==1)+I(pclass==2)+
I(sex=="female")+
age+
fare+
I(embarked=="C")+I(embarked=="Q")+
I(Title=="Master") + I(Title=="Miss") + I(Title=="Mr") +
I(family_size_discreto == "solteiro")+I(family_size_discreto == "grande")+
I(child=="crianca")+
I(Mother== "mae"),
data = data_train2, family=binomial(link = "probit"))
cast2 <- predict(Probit, newdata = data_test2, type = "response")
prediction2 <- table(data_test2$survived, cast2 > 0.5)
TP2 <- prediction2["1", "TRUE"]
FP2 <- prediction2["0", "TRUE"]
TN2 <- prediction2["0", "FALSE"]
FN2 <- prediction2["1", "FALSE"]
accuracy2[i] = (TP2 + TN2)/(TP2 + TN2 + FP2 + FN2)}
paste("A Acurácia média em", toString(k_folds), "folds é", toString(round(mean(accuracy2), 4)))## [1] "A Acurácia média em 5 folds é 0.8197"Abaixo, segue o resumo do modelo Probit:
summary(Probit)##
## Call:
## glm(formula = survived ~ I(pclass == 1) + I(pclass == 2) + I(sex ==
## "female") + age + fare + I(embarked == "C") + I(embarked ==
## "Q") + I(Title == "Master") + I(Title == "Miss") + I(Title ==
## "Mr") + I(family_size_discreto == "solteiro") + I(family_size_discreto ==
## "grande") + I(child == "crianca") + I(Mother == "mae"), family = binomial(link = "probit"),
## data = data_train2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3264 -0.6162 -0.4117 0.5882 3.3519
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) -9.756e-01 1.904e-01 -5.125
## I(pclass == 1)TRUE 1.107e+00 1.471e-01 7.528
## I(pclass == 2)TRUE 4.968e-01 1.265e-01 3.928
## I(sex == "female")TRUE 1.923e+00 1.869e-01 10.286
## age -1.705e-02 4.970e-03 -3.431
## fare 1.561e-07 1.339e-07 1.166
## I(embarked == "C")TRUE 2.847e-01 1.292e-01 2.204
## I(embarked == "Q")TRUE 4.010e-02 1.731e-01 0.232
## I(Title == "Master")TRUE 1.447e+00 2.767e-01 5.230
## I(Title == "Miss")TRUE -3.483e-01 2.113e-01 -1.648
## I(Title == "Mr")TRUE NA NA NA
## I(family_size_discreto == "solteiro")TRUE 7.346e-02 1.196e-01 0.614
## I(family_size_discreto == "grande")TRUE -1.218e+00 2.305e-01 -5.282
## I(child == "crianca")TRUE -2.093e-01 2.004e-01 -1.044
## I(Mother == "mae")TRUE 4.461e-01 2.804e-01 1.591
## Pr(>|z|)
## (Intercept) 2.98e-07 ***
## I(pclass == 1)TRUE 5.14e-14 ***
## I(pclass == 2)TRUE 8.55e-05 ***
## I(sex == "female")TRUE < 2e-16 ***
## age 0.000602 ***
## fare 0.243759
## I(embarked == "C")TRUE 0.027531 *
## I(embarked == "Q")TRUE 0.816791
## I(Title == "Master")TRUE 1.69e-07 ***
## I(Title == "Miss")TRUE 0.099250 .
## I(Title == "Mr")TRUE NA
## I(family_size_discreto == "solteiro")TRUE 0.539207
## I(family_size_discreto == "grande")TRUE 1.28e-07 ***
## I(child == "crianca")TRUE 0.296367
## I(Mother == "mae")TRUE 0.111594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1392.63 on 1046 degrees of freedom
## Residual deviance: 891.23 on 1033 degrees of freedom
## AIC: 919.23
##
## Number of Fisher Scoring iterations: 54.5 Conclusão: Comparando os modelos:
Agora faremos a comparação final dos três modelos:
###COMPARANDO OS TRÊS MODELOS
RESULT_LOGIT <- paste("Acurácia média", toString(k_folds), "folds", toString(round(mean(accuracy), 4)))
RESULT_PROBIT <- paste("Acurácia média", toString(k_folds), "folds", toString(round(mean(accuracy2), 4)))
RESULT_ARVORE <-paste("Acurácia média", toString(k_folds), "folds", toString(round(mean(accuracy_tree2), 4)))RESULT_LOGIT ## [1] "Acurácia média 5 folds 0.8182"RESULT_PROBIT ## [1] "Acurácia média 5 folds 0.8197"RESULT_ARVORE## [1] "Acurácia média 5 folds 0.8143"A conclusão é que a acurácia dos modelos é muito parecida, não havendo diferença significtiva que justifique a escolha de um deles.