Escopo

Quem gosta de conhecer bem os alimentos que consome, de manter dietas e um estilo de vida saudável já deve ter ouvido falar nos grupos alimentares. Sua divisão ocorre por categoria dos alimentos e diverge um pouco entre estudiosos da área, enquanto alguns dizem que são 8 grupos, outros afirmam ser 6, ou até mesmo 5. O que importa de fato são os fatores que compõem cada alimento Afinal, são as características calóricas, proteicas, lipídicas, entre outras que definem se o alimento deve ou não entrar na dieta.

Posto isso, me deparo com um case interessante. Minha irmã foi ao nutricionista, e este te passou uma lista contendo 20 alimentos. Junto à lista veio a indicação: “Consuma apenas um único alimento de cada grupo por dia”. Chegando em casa minha irmã se deu conta que a literatura sobre os grupos eram um pouco contraditória e resolveu ligar para o nutricionista. Contudo, como todo bom e velho cientista de dados curioso, resolvi me adiantar à resposta do nutricionista e ajudar minha irmã a dividir os alimentos em grupos.

O resultado ficou tão legal que resolvi expandir a modelagem para o caso de minha irmã querer substituir o alimento do lista por algum outro, fora da lista, mas que pertencesse ao mesmo grupo. Dá uma conferida nesse passo a passo, está bem didático e é um bom norte para quem pensa em trabalhar na área de ciência de dados.

Análise dos Dados

Aqui, mostrarei um breve descritivo dos dados, tentando ao máximo separá-lo em grupos.

library(foreign)
library(skimr)
library(haven)
library(stringr)
library(readr)
dados = read_delim("alimentos.csv", ";", escape_double = FALSE, trim_ws = TRUE)
colnames(dados) = c('Alimento','Energia_kcal','Proteinas_g','Lipidos_g','Calcio_mg','Ferro_mg')
dados = as.data.frame(dados)

A lista de minha irmã está descrita na variável Alimento. As demais variáveis dizem respeito às informações que busquei a cerca de cada um deles.

kable(dados, format = "markdown", digits = 2)
Alimento Energia_kcal Proteinas_g Lipidos_g Calcio_mg Ferro_mg
Azeite 900 0.0 100.0 0.1 0.05
Manteiga 770 0.0 85.0 13.0 0.20
Pescada 85 19.0 1.0 25.0 0.90
Vaca 208 18.0 15.0 12.0 1.50
Frango 158 20.0 8.5 18.0 1.80
Leite 57 3.0 3.0 126.0 0.10
Iogurte 59 3.2 3.2 125.0 0.20
Q. Flamengo 316 26.0 23.2 800.0 0.80
Q. Serra 392 26.0 32.0 800.0 1.20
Arroz 350 7.5 0.5 10.0 0.50
Pao 258 7.0 0.6 24.0 1.60
Feijao 290 20.0 1.2 170.0 6.50
Acucar 400 0.0 0.0 15.0 1.00
Massas 365 10.0 0.5 20.0 1.00
Alface 22 1.8 0.2 70.0 1.50
Cebola 22 0.9 0.2 31.0 0.50
Espinafres 22 2.6 0.9 104.0 3.60
Cenoura 22 0.6 0.0 104.0 3.60
Batata 90 2.5 0.0 9.0 0.20
Couve 30 2.9 0.5 234.0 1.80

Primeiramente, utilizei a matriz de correlação para verificar uma possível correlação linear entre as variáveis. Quando mais escuras nas tonalidades vermelho ou azul, maior a magnitude da relação linear entre as variáveis. Bem como, quanto mais inteiro o gráfico de pizza, maior a correlação linear.

Como o interesse não é reduzir a dimensionalidade, mas sim separar grupos, quanto menos correlacionadas estiverem as variáveis mais facilitado o trabalho de dividi-las em grupos. Contudo, um comportamento linearmente similar similar entre entre duas ou três variáveis pode ajudar a reduzir as cinco covariáveis à apenas duas. O que facilitaria no momento de análise visual acerca das fronteiras de decisão.

library(corrplot)
corrplot(cor(dados[,-1]), type="upper", method="pie")

Turning do Número de Clusters

Feito isso, a primeira ideia para definição do número de grupos foi um turning da performance dos clusters serapados por K-Means. A ideia foi definir um número de cluster tal que a Soma do Quadrado dos Resíduos (SQR) fose mínima, mas lembrando sempre do Tradeoff Vício-Variabilidade, afinal, não queremos overffiting no nosso modelo.

Posto isso, perceba pelo gráfico abaixo, que a partir de 5 clusters, o ganho no decréscimo da SQR foi mínimo. Posto isso, vamos supor para os próximos passos o número ideal de 5 clusters.

x = cbind("Qt. Energetica (kcal)" = dados[,2], "Qt. de Proteinas (g)" = dados[,3],
          "Qt. de Lipidios (g)" = dados[,4], "Qt. de Calcio (mg)" = dados[,5],
          "Qt. de Ferro (mg)" = dados[,6])
chute.num.clusters = function(n,x){
contador = c(rep(NA,n)) 
for (i in 1:n) {
  saida = kmeans(x, centers = i, nstart = 20)
  contador[i] <- saida$tot.withinss
}
modelo = plot(1:n, contador, type = "b", main = "Turning do Numero de Clusters",col = 'tomato',
     xlab = "Numero de Clusters", 
     ylab = "Soma dos Quadrados dos Residuos por Grupos")
}
chute.num.clusters(n=10,x=x)

Criação de Nova Variável por Grupo a priori

Posto isso, criei 5 grupos, dado informações a priori sobre a classificação de alimentos, foram eles: proteina, graos, vegetais, laticinios e acucar.

proteina = c("Pescada","Vaca","Frango")
graos = c("Arroz","Feijao")
vegetais = c("Alface","Cebola","Espinafres","Cenoura","Batata","Couve")
laticinios= c("Leite","Iogurte","Q. Flamengo","Q. Serra","Pao","Massas","Manteiga","Azeite")
acucar = c("Acucar")
posicao.preteinas = match(proteina,dados[,1])
posicao.vegetais = match(vegetais,dados[,1])
posicao.graos = match(graos,dados[,1])
posicao.laticinios = match(laticinios,dados[,1])
posicao.acucar = match(acucar,dados[,1])
dados$Tipo = c(rep(NA,nrow(dados)))
dados$Tipo[posicao.preteinas] = rep('proteina', length(posicao.preteinas))
dados$Tipo[posicao.graos] = rep('grao', length(posicao.graos))
dados$Tipo[posicao.vegetais] = rep('vegetal', length(posicao.vegetais))
dados$Tipo[posicao.laticinios] = rep('laticinio', length(posicao.laticinios))
dados$Tipo[posicao.acucar] = rep('acucar', length(posicao.acucar))
dados$Tipo = as.factor(dados$Tipo)

Ao fazer uma breve análise exploratória univariada dos dados, a partir da nova variável criada com base nos grupos, espera-se observar diferenças substanciais que tornem coerente essa divisão.

boxplot(dados[,2]~dados$Tipo, ylim= c(0,900), main = 'Grupo Alimentar vs Qtd. Energia (kcal)', xlab= 'Grupo',ylab = '', col=c('lightblue','tomato','orange', 'green'))

boxplot(dados[,3]~dados$Tipo, ylim= c(0,30), main = 'Grupo Alimentar vs Qtd. Proteina (g)', xlab= 'Grupo',ylab = '', col=c('lightblue','tomato','orange', 'green'))

boxplot(dados[,4]~dados$Tipo, ylim= c(0,100), main = 'Grupo Alimentar vs Qtd. Lipidio (g)', xlab= 'Grupo',ylab = '', col=c('lightblue','tomato','orange', 'green'))

boxplot(dados[,5]~dados$Tipo, ylim= c(0,800), main = 'Grupo Alimentar vs Qtd. Calcio (mg)', xlab= 'Grupo',ylab = '', col=c('lightblue','tomato','orange', 'green'))

boxplot(dados[,6]~dados$Tipo, ylim= c(0,7), main = 'Grupo Alimentar vs Qtd. Ferro (mg)', xlab= 'Grupo',ylab = '', col=c('lightblue','tomato','orange', 'green'))

Percebe-se uma certa diferença comportamental entre os grupos, o que enfatiza a coerência na divisão. Contudo, não me limitarei a ela.

Modelagem Via Clusterização

Aqui, utilizarei tanto a informação a priori acerca da divisão dos grupos, quanto deixarei a cargo do próprio modelo a forma que será dada essa divisão, fixando apenas o número de clusters definido na etapa de tunagem.

library(ggfortify)
library(lfda)
library(stats)
library(cluster)
xis = dados[,-c(1,7)]
pca = autoplot(prcomp(x), data = dados, frame = TRUE,colour = "Tipo", main = "Metodo PCA")
kmeans = autoplot(kmeans(xis, 5), data = dados, label.size = 5, frame = TRUE,
                  colour = "Tipo", main = " Metodo K-Means")
clara = autoplot(clara(xis, 5), frame = TRUE, main = "Metodo Clara")
fanny = autoplot(fanny(xis, 5), frame = TRUE, main = "Metodo Fanny")
pam = autoplot(pam(xis, 5), frame = TRUE,  main = "Metodo Pam")

Modelos com Consideração a Priori

Aqui, utilizei a variável criada a partir dos grupos definidos a partir de pesquisa conceitual para modelar os dados por dois métodos: PCA e K-Means. PCA, ou Componentes principais é o nome dado ao conjunto de valores de variáveis linearmente não correlacionadas, convertidas a partir de observações de variáveis linearmente correlacionadas. Já, o K-Means é um método de Clusterização que objetiva particionar n observações dentre k grupos onde cada observação pertence ao grupo mais próximo da média. Dividindo o espaço em um diaframa de Voroni, o K-Means é melhor quanto maior a diferença entre as médias de cada classe.

pca

kmeans

Perceba que a utilização da informação a priori não resultou em uma divisão ‘ótima’, uma vez que a fronteira de decisão dos laticíneos engloba açúcares e lipídios. O que pode refletir em erros de previsão para novos dados.

Modelos sem Consideração da Informação a Priori

Imaginando que os pontos não identificados pelos métodos PCA e K-Means poderiam ser possíveis outliers, testei um algoritmo mais robusto à eles, em comparação com o K-Means, porque minimiza uma soma de “dessemelhanças” aos pares, em vez de uma soma de distâncias euclidianas quadradas, este método é o PAM, ou K-Medoids. Este, pode ser definido como o objeto de um cluster cuja diferença média para todos os objetos no cluster é mínima, ou seja, é o ponto mais centralmente localizado no cluster. Ou seja, diferentemente do K-Means, o centro de um cluster não é necessariamente um dos pontos de dados de entrada (é a média entre os pontos no cluster).

pam

Percebe-se então, que o PAM apresentou ganho significativo na clusterização, quando comparado aos métodos anteriores, porém deixou muita áreas difusas para previsões futuras. Por isso tentei o Método CLARA. Este, funciona agrupando uma amostra do conjunto de dados e atribui todos os objetos no conjunto de dados a esses clusters.

clara

Contudo, o método CLARA não apresentou ganho quando comparado ou PAM. Por isso, resolvi tester um último método. o Método Fanny. Este, é um algoritmo de agrupamento suave ou difuso, em que cada nó no gráfico está associado a um coeficiente de associação, indicando o grau de pertencimento de cada nó a diferentes clusters. A idéia então é expandir a região de decisão dos grupos, diminuindo o máximo possível de regiões difusas.

fanny

O que observou-se foi que de fato, o método Fanny, apresentou uma melhora, ainda que sensível, na disposição das suas fronteiras de decisão.

Definido o modelo de clusterização, estimamos a forma em que cada alimento foi categorizado conforme o modelo selecionado. Feito isso, criei um novo banco de dados, com os alimentos já agrupados.

modelo = fanny(xis, 5)
agregados = scale(xis)
novos.dados =  data.frame(agregados, modelo$clustering)
library(mclust)
fit = Mclust(novos.dados)
GRUPO.PREDITO= as.factor(predict(fit)$classification)
dados.agrupados= dados[,-c(1,7)]
dados.agrupados$GRUPO.PREDITO = GRUPO.PREDITO
kable(dados.agrupados, format = "markdown", digits = 2)
Energia_kcal Proteinas_g Lipidos_g Calcio_mg Ferro_mg GRUPO.PREDITO
900 0.0 100.0 0.1 0.05 1
770 0.0 85.0 13.0 0.20 1
85 19.0 1.0 25.0 0.90 2
208 18.0 15.0 12.0 1.50 2
158 20.0 8.5 18.0 1.80 2
57 3.0 3.0 126.0 0.10 3
59 3.2 3.2 125.0 0.20 3
316 26.0 23.2 800.0 0.80 1
392 26.0 32.0 800.0 1.20 1
350 7.5 0.5 10.0 0.50 4
258 7.0 0.6 24.0 1.60 4
290 20.0 1.2 170.0 6.50 1
400 0.0 0.0 15.0 1.00 4
365 10.0 0.5 20.0 1.00 4
22 1.8 0.2 70.0 1.50 3
22 0.9 0.2 31.0 0.50 3
22 2.6 0.9 104.0 3.60 1
22 0.6 0.0 104.0 3.60 1
90 2.5 0.0 9.0 0.20 3
30 2.9 0.5 234.0 1.80 3

Modelos para Previsão de Novos Alimentos

Com esse banco de dados, a ideia agora é utilizar modelos supervisionados possibilitando prever a categoria dos alimentos, conforme agrupamento realizado na fase anterior. Para isso, utilizarei o índice kappa como medida de performance para seleção do modelo ótimo.

dados = dados.agrupados
varResp = 'GRUPO.PREDITO'
Medidas = function(mc){
  PE = (sum(mc[1,])*sum(mc[,1])+sum(mc[2,])*sum(mc[,2])+
    sum(mc[3,])*sum(mc[,3])+sum(mc[4,])*sum(mc[,4]))/((sum(mc[1,])+sum(mc[2,])+sum(mc[3,])+sum(mc[4,]))^2)
  PA = sum(diag(mc))/(sum(mc[1,])+sum(mc[2,])+sum(mc[3,])+sum(mc[4,]))
  Kappa = (PA-PE)/(1-PE)
  return(list(Kappa = Kappa))
}
formula = as.formula(paste(varResp, '~ .'))

FormFac = function(formula){
  aux = as.character(formula)[2]
  return(as.formula(paste('factor(',aux, ')', '~ .')))
}
AuxSaida = function(pred, teste, formula, alg){
  ind = which(colnames(teste) == as.character(formula)[2])
  mc = matrix(NA,4,4)
  mc[1,1] = sum(teste[, ind]==1 & pred ==1)
  mc[1,2] = sum(teste[, ind]==2 & pred ==1)
  mc[1,3] = sum(teste[, ind]==3 & pred ==1)
  mc[1,4] = sum(teste[, ind]==4 & pred ==1)
  mc[2,1] = sum(teste[, ind]==1 & pred ==2)
  mc[2,2] = sum(teste[, ind]==2 & pred ==2)
  mc[2,3] = sum(teste[, ind]==3 & pred ==2)
  mc[2,4] = sum(teste[, ind]==4 & pred ==2)
  mc[3,1] = sum(teste[, ind]==1 & pred ==3)
  mc[3,2] = sum(teste[, ind]==2 & pred ==3)
  mc[3,3] = sum(teste[, ind]==3 & pred ==3)
  mc[3,4] = sum(teste[, ind]==4 & pred ==3)
  mc[4,1] = sum(teste[, ind]==1 & pred ==4)
  mc[4,2] = sum(teste[, ind]==2 & pred ==4)
  mc[4,3] = sum(teste[, ind]==3 & pred ==4)
  mc[4,4] = sum(teste[, ind]==4 & pred ==4)
  rs = Medidas(mc)
  rs[['Classificador']] = alg    
  return(do.call(data.frame, rs))  
}

AuxSaida.bayes = function(pred, teste, formula, alg){
  ind = which(colnames(teste) == as.character(formula)[2])
  mc = matrix(NA,4,4)
  mc[1,1] = sum(teste[, ind]==1 & pred ==1)
  mc[1,2] = sum(teste[, ind]==2 & pred ==1)
  mc[1,3] = sum(teste[, ind]==3 & pred ==1)
  mc[1,4] = sum(teste[, ind]==4 & pred ==1)
  mc[2,1] = sum(teste[, ind]==1 & pred ==2)
  mc[2,2] = sum(teste[, ind]==2 & pred ==2)
  mc[2,3] = sum(teste[, ind]==3 & pred ==2)
  mc[2,4] = sum(teste[, ind]==4 & pred ==2)
  mc[3,1] = sum(teste[, ind]==1 & pred ==3)
  mc[3,2] = sum(teste[, ind]==2 & pred ==3)
  mc[3,3] = sum(teste[, ind]==3 & pred ==3)
  mc[3,4] = sum(teste[, ind]==4 & pred ==3)
  mc[4,1] = sum(teste[, ind]==1 & pred ==4)
  mc[4,2] = sum(teste[, ind]==2 & pred ==4)
  mc[4,3] = sum(teste[, ind]==3 & pred ==4)
  mc[4,4] = sum(teste[, ind]==4 & pred ==4)
  rs = Medidas(mc)
  rs[['Classificador']] = alg    
  return(do.call(data.frame, rs))  
}

Funções de Data Splitting

Aqui, apresento funções referentes ao processo de data splitting. Neste trabalho, fiz a validação baseando-se na validação cruzada k-fold repetida. Para um k=10, e 10 repetições.

Contudo, assim como na seção anterior, trago outras alternativas a quem for de interesse a replicar esse trabalho. Desta forma, além da própria validação k-fold, sem repetições, uma outra possibilidade é o uso do Holdout, implementado através da função denominada Holdout, em que o parâmetro p, indica o percentual utilizado na amostra de treino.

###### DATA SPLITTING 
Holdout = function(dados, p = 0.7){
  n_treino = ceiling(dim(dados)[1]*p)
  ind = c(rep('treino', n_treino), rep('teste', dim(dados)[1] - n_treino) )
  ind = sample(ind)
  treino = dados[ind == 'treino', ]
  teste  = dados[ind == 'teste' , ]
  return(list(treino = treino, teste = teste))
}


ValidCruzada = function(Bayesiano = 'n', cls, dados, formula, kfolds = 10, ...){
  library(dismo)
  library(bnlearn)
  dados.discretizados = discretize(dados,breaks = 3)
  if (Bayesiano == 's') {dados = dados.discretizados}
  id = kfold(1:nrow(dados), kfolds)
  vcs = data.frame()
  for(i in 1:kfolds){
    treino = dados[id != i, ]
    teste = dados[id == i, ]
    kcls = cls(treino, teste, formula,...) 
    vcs = rbind(vcs, kcls)
  }
  vcs = vcs[, c(ncol(vcs), 1:(ncol(vcs) - 1))]
  avg = aggregate(. ~ Classificador, data = vcs, FUN = mean)
  std = aggregate(. ~ Classificador, data = vcs, FUN = sd)
  return(list(Media = avg, Desvio = avg, Modelos = vcs))
}
ValCruzRep = function(Bayesiano = 'n',cls, dados, formula, kfolds = 10, reps = 10,...){
  x = data.frame()
  for(i in 1:reps) x = rbind(x, ValidCruzada(Bayesiano = Bayesiano,cls, dados, formula, kfolds,...)$Media) 
  avg = aggregate(. ~ Classificador, data = x, FUN = mean)
  std = aggregate(. ~ Classificador, data = x, FUN = sd)
  return(list(Media = avg, Desvio = std, Modelos = x))
}

Funções dos Modelos

Nesta seção, apresento as funções de cada um dos modelos estudados: Regressão Logística, LDA, Floresta Aleatória, Naive Bayes, KDB e SVM. Em que, para os modelos Floresta Aleatória, KDB E SVM, apresento as funções acerca do processo de turning sobre seus parâmetros: o número de árvores, o número de dependências entre as covariáveis explicativas e o kernel, respectivamente.

########## Logistica 
RegLog = function(treino, teste, formula){ 
  cls = glm(formula, data = treino, family = 'binomial')
  pred = predict(cls, newdata = teste, type = 'response') 
  pred = ifelse(pred > 0.5, 1, 0) 
  pred = as.numeric(pred)
  pred = as.factor(pred)
  return(pred)
}

medRegLog = function(treino, teste, formula){
  alg = 'Regressao Logistica' 
  pred = RegLog(treino, teste, formula)
  return(AuxSaida(pred, teste, formula, alg))
}


########## LDA
LDA = function(treino, teste, formula){
  library(MASS) 
  cls = lda(formula, data = treino)
  pred = predict(cls, newdata = teste)$class 
  pred = as.numeric(pred) - 1
  return(pred)
}

medLDA = function(treino, teste, formula){ 
  alg = 'Analise Discriminante Linear' 
  pred = LDA(treino, teste, formula)
  return(AuxSaida(pred, teste, formula, alg))
}


########## Naive Bayes 
NB = function(treino, teste,varResp, formula){
  library(bnclassify)
  cls = bnc('nb', varResp, treino, smooth = 1)
  pred = predict(cls, teste)
  return(pred)
}

medNB = function(treino, teste,formula){ 
  alg = 'Naive Bayes' 
  pred = NB(treino, teste,varResp, formula)
  return(AuxSaida.bayes(pred, teste, formula, alg))
}


########## KDB
KDB = function(treino, teste,varResp, formula,K=2){
  library(bnclassify)
  cls = kdb(varResp,treino, k=10, kdbk = K, epsilon = 0.01, smooth = 0,cache_reset = NULL)
  cls  = lp(cls, treino, smooth = 1)
  pred = predict(cls,teste,prob = FALSE)
  return(pred)
}

medKDB = function(treino, teste,formula,K=2){ 
  alg = paste('KDB')
  pred = KDB(treino, teste,varResp, formula,K=K)
  return(AuxSaida.bayes(pred, teste, formula, alg))
}

 
# Turning de K dependentes do KDB, via Validação Cruzada K=10 folds
KDBmelhorK = function(dados, formula, inicio = 1, fim = (dim(dados)[2]-2), metrica = 'Kappa', ...){
  knns = data.frame()
  s = seq(inicio, fim, by = 1)
  for(k in s) knns = rbind(knns, ValidCruzada(Bayesiano = 's',cls = medKDB,dados, formula, kfolds = 10,K=k)$Media)
  ordem = order(knns[metrica])
  melhorK = s[tail(ordem, 1)]
  ord = knns[ordem, ]
  plot(s, knns[[metrica]], pch = 17, xlab = 'K- Dependentes', main = paste('Melhor K - metrica:', metrica),
       ylab = metrica, xaxt = "n", col = 'blue') 
  axis(1, at = s)
  segments(x0 = melhorK, y0 = 0, x1 = melhorK, y1 = as.numeric(tail(ord[metrica], 1)), 
           lty = 'dotted', col = 'gray50', lwd = 2)
  return(melhorK)
}

SVM = function(treino, teste, formula,KERNEL){
  library(kernlab)
  cls = ksvm(formula, data = data.frame(treino),type='C-svc', kernel= KERNEL)
  pred = predict(cls, newdata = data.frame(teste))
  return(pred)
}
medSVM = function(treino, teste, formula,KERNEL){
  alg = paste0('SVM-',KERNEL)
  pred = SVM(treino, teste, formula,KERNEL)
  return(AuxSaida(pred, teste, formula, alg))
}

# Turning de Kernels, via Validação Cruzada K=10 folds
SVMmelhorKernel = function(dados, formula, metrica = 'Kappa', ...){
  svmturning = data.frame()
  s = c("rbfdot","polydot","vanilladot")
  for(i in 1:length(s)){
    KERNEL = s[i]
    svmturning = rbind(svmturning, ValidCruzada(Bayesiano = 'n',medSVM, dados, formula, KERNEL = KERNEL)$Media)
  }
  return(kable(svmturning))
}

Floresta = function(treino, teste, formula,ARVORES){
  library(randomForest) 
  cls = randomForest(formula, data = data.frame(treino),importance= TRUE,ntree= ARVORES,
                     proximity=TRUE)
  pred = predict(cls, newdata = data.frame(teste), type = "class")
  return(pred)
}
medFloresta = function(treino, teste, formula,ARVORES){ 
  alg = 'Floresta Aleatoria' 
  pred = Floresta(data.frame(treino), data.frame(teste), formula,ARVORES)
  return(AuxSaida(pred, teste, formula, alg))
}

# Turning de Arvores, via Validação Cruzada K=10 folds
RFmelhorARVORE = function(dados, formula, inicio = 1, fim = (dim(dados)[1]/2), metrica = 'Kappa', ...){
  RFturning = data.frame()
  s = seq(inicio, fim, by = 1)
 for(k in s) RFturning = rbind(RFturning, ValidCruzada(Bayesiano = 'n',medFloresta, dados, formula, ARVORES = k, ...)$Media)

  ordem = order(RFturning[metrica])
  melhorARVORE = s[tail(ordem, 1)]
  ord = RFturning[ordem, ]

  plot(1:length(s), RFturning[[metrica]], pch = 17, xlab = 'Arvores', main = paste('Melhor Arvore - metrica:', metrica),
       ylab = metrica, xaxt = "n", col ='blue') 
  axis(1, at = s)
  segments(x0 = melhorARVORE, y0 = 0, x1 =melhorARVORE, y1 = as.numeric(tail(ord[metrica], 1)), 
           lty = 'dotted', col = 'gray50', lwd = 2)
  
  return(melhorARVORE)
}

Turning dos KDB, Floresta e SVM

Aqui, apliquei as funções de tunagem apresentadas acima, a fim de escolher o melhor K para o modelo KDB, variando o número de k-dependentes considerados, bem como o melhor número de Árvores e o melhor Kernel, para os modelos Floresta Aleatória e SVM, respectivamente. Lembre que esses parâmetros representam a complexidades desses modelos, e portanto merece atenção quanto ao tradeoff viés-variância. Isto evita possíveis situações de over(under)fitting.

#TURNING KDB 
KDBmelhorK(dados, formula, 1 ,(dim(dados)[2]-2),'Kappa', kfolds = 10)

## [1] 4
  SVMmelhorKernel(dados, formula, metrica = 'Kappa')
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters
Classificador Kappa
SVM-rbfdot 1.0000000
SVM-polydot 1.0000000
SVM-vanilladot 0.8333333 
RFmelhorARVORE(dados, formula, inicio = 1, fim = (dim(dados)[1]/2), metrica = 'Kappa')
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin

## [1] 9

Desempenho dos Modelos

Por fim, aplicando a função ValCruzRep, ilustrada anteriormente, construímos, com base no índice Kappa, a tabela de desempenho dos modelos estudados. Observe que o desempenho em todos os modelos foi muito bom, o que evidencia uma clusterização adequada dos alimentos.

Além disso, ao analisar a tabela abaixo, destaco que, entre os modelos estudados, o de melhor performance preditiva foi o Modelo SVM com Kernel Gaussiano. Vale a ressalva interpretativa que os modelos para índice Kappa menores que 0.25, apresentam resultado piores que a aleatoriedade. Uma vez que, para 4 classes de predição, um modelo naive acertaria 25% dos casos.

tabela.desempenho = matrix(NA,6,2)
colnames(tabela.desempenho) = c('Modelo', 'Kappa')
saida = ValCruzRep(Bayesiano = 'n',cls = medRegLog,dados, formula, kfolds = 10)
tabela.desempenho[1,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medLDA,dados, formula, kfolds = 10)
tabela.desempenho[2,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 's',cls = medNB,dados, formula, kfolds = 10)
tabela.desempenho[5,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 's',cls = medKDB,dados, formula, kfolds = 10, K= 1)
tabela.desempenho[6,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medSVM,dados, formula, kfolds = 10, KERNEL='rbfdot')
tabela.desempenho[3,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medFloresta,dados, formula, kfolds = 10, ARVORES= 10)
tabela.desempenho[4,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
tabela.desempenho= tabela.desempenho[order(tabela.desempenho[,2]),]


library('formattable')
tabela.desempenho = data.frame(tabela.desempenho)
formattable(tabela.desempenho, list(Acuracia =color_tile("orange", "red")))
Modelo Kappa
Analise Discriminante Linear -0.085
Regressao Logistica 0
KDB 0.199
Naive Bayes 0.262
Floresta Aleatoria 0.564
SVM-rbfdot 0.928