Qual a melhor partição dos meus dados, ou seja, qual o melhor número de grupos que meus dados sugerem?

• hennig2013.pdf - Leitura importante sobre suposições em classificações não supervisionadas (autor do pacote fpc do R)

A estratégia aqui implementada é bem simples:

1. Defina uma sequência de números de grupos (k) a ser testada
2. Para cada valor de (k):
   1. define agrupamentos para k por vários métodos diferentes (kmeans, pam, specc, mclust) - com dados morfológicos
   2. testa cada agrupamento com svm,bayes e lda (pode usar os mesmos dados que gerou os grupos ou outro conjunto de dados, e.g. NIR)
   3. faz uma figura mostrando os resultados desses testes (para cada K e metodo mostra a proporcao de amostras preditas corretamente)
3. O número ótimo de K é o valor máximo a partir do qual a taxa de acerto máxima diminui (ponto de inflexão ou diretamente pelo critério aceitável, e.g. 95% das amostras). Ver tomada de decisão flexível abaixo.
4. O agrupamento final, segundo o K ótimo, é refeito e testada e as amostras são categorizadas em 3 classes:
   1. GRUPO CORE - amostras preditas corretamente com prob>0.95 (classe 1)
   2. INCERTEZAS (zona de tolerância entre grupos):
      1. amostras preditas corretamente mas com prob⇐0.95 (classe 0)
   3. SEM CLASSE
      1. amostras reclassificadas num grupo diferente do atribuído inicialmente
      2. amostras que em métodos igualmente bons são classificadas diferentes
5. O resultado é uma tabela com o(s) agrupamento(s), as probabilidades de pertenecimento em cada grupo e a(s) classificacao(oes) acima. Isso é salvo num arquivo que tem no nome classificacaoAposterioriK.

• funcoesnecessarias.r - precisa dos pacotes e funções chamados neste arquivo
• dadosmorfonmdsvarsel.csv - arquivo de dados usado no exemplo
• otimongruposkmeanstest.pdf - exemplo da figura gerada

rm(list=ls())
library(mclust)
 
#funcoes e pacotes necessarios
source("funcoesNecessarias.r")
 
#para salvar resultados
pathparadados = "dadosParaAnalises"
pathparafiguras = "figuras"
pathparatabelas = "tabelas"
#dir.create(pathparadados)
#dir.create(pathparafiguras)
#dir.create(pathparatabelas)
 
#pega os dados
#o exemplo testado é uma matriz de NMDS gerada a partir
#de distancias morfológicas
arq = "dadosParaAnalises/dadosMORFOnmdsvarsel.csv"
dados.morfo = read.table(file=arq,sep='\t',as.is=F,row.names = 1,header=T)
 
#PASSO 01 - QUAL O NÚMERO DE GRUPOS?
#qual a porcentagem de acerto usando testa.classificacao
#para cada k e metodo de particao (mclust,kmeans, pam ou specc)?
#ver help da funcao clusterboot() que e uma alternativa ao
#implementado abaixo
#a lógica neste neste é que o melhor k (e método de particao)
#que indique agrupamentos das amostras
#é aquele em que há maior taxa de acerto numa validacao cruzada
 
#a validacao pode ser feita com os mesmos dados que geram os grupos
#ou com dados diferentes
dado.grupos = dados.morfo
dado.teste = dados.morfo #pode usar NIR aqui por exempo, se tiver
#ambos esses dados precisam ter o mesmo número
#de linhas e as amostras na mesma ordem
#garante isso
dado.teste = dado.teste[row.names(dado.grupos),]
 
#qual a sequencia de valores de K a serem testados (k=número de grupos)
kseq = seq(16,2,by=-1) 
#para cada valor de K, gera diferentes agrupamentos para diferentes metodos
ltotest = NULL
for(oka in kseq) {
  cat("0\14")
  cat(paste("Gerando K=",oka,"grupos"))
  tkn = speccCBI(dado.grupos,k=oka)
  ltotest[["speccCBI"]][[as.character(oka)]] = tkn$partition
  clb1 = Mclust(dado.grupos,G=oka)
  ltotest[["Mclust"]][[as.character(oka)]] = clb1$classification
 
  clb1 = pam(dado.grupos,k=oka,cluster.only=T)
  ltotest[["pam"]][[as.character(oka)]] = clb1
 
  clb1 = kmeans(dado.grupos, center=oka,iter.max=10000,nstart=1000)
  ltotest[["kmeans"]][[as.character(oka)]] = clb1$cluster
}
 
str(ltotest[[2]])
#agora testa a estabilidade dos agrupamentos
#isso pode demorar bastante - para cada grupo acima ira testar a classificacao dos dados 
bestks.tb = data.frame(FUN=NA,K=NA,MAX.TAXA.ACERTO=NA,SVM=NA,BAYES=NA,LDA=NA,SVM.ALL=NA,BAYES.ALL=NA,LDA.ALL=NA)
bestks = NULL
for(l in 1:length(ltotest)) {
  lt = ltotest[[l]]
  nlt = names(ltotest)[l]
  for(n in 1:length(lt)) {
    print(paste(names(ltotest)[l],"k=",names(ltotest[[l]])[n]))
    #pega o agrupamento
    cl = ltotest[[l]][[n]]
    oka = names(ltotest[[l]])[n]
    #copia o dado
    dd = dado.teste
    #adiciona o agrupamento gerado
    rn = row.names(dd)
    dd$GROUPING = as.factor(paste("GRUPO",cl))
 
    #testa o agrupamento com lda loo
    oteste = testa.classificacao(dad=dd,colgrp="GROUPING",analises=c("svm","bayes","lda"),method=c("loo","perc")[1],perc=15)
    #pega o fit do teste (maior taxa de acerto total)
    br =oteste$BRUTOS
 
    #pega melhor predicao
    vl = br$MELHOR.PREDICAO.PROB>0.95 & br$ACERTO=='acerto' & !is.na(br$ACERTO)
    ovv = table(as.vector(br[vl,"ANALISE"]))
    ovv1 = ovv/length(unique(br$ID))
    cll = c("BAYES","LDA","SVM")
    mcll = cll[!cll%in%names(ovv1)]
    rz = c(FUN=nlt,K=oka,MAX.TAXA.ACERTO=max(ovv1),ovv1)
    if (length(mcll)>0) {rz[mcll] = NA}
 
    #pega acertos totais
    vl = br$ACERTO=='acerto'
    ovv = table(as.vector(br[vl,"ANALISE"]))
    ovv = ovv/length(unique(br$ID))
    cll = c("BAYES","LDA","SVM")
    mcll = cll[!cll%in%names(ovv)]
    if (length(mcll)>0) {ovv[mcll] = NA}
    ovv = ovv[cll]
    names(ovv) = c("BAYES.ALL","LDA.ALL","SVM.ALL")
    rz = c(rz,ovv)
    rz = rz[colnames(bestks.tb)]
    print(paste(round(max(ovv,na.rm=T),4)," taxa de acerto máxima para K=",oka))
    #salva o resultado do teste para k
    bestks.tb = rbind(bestks.tb,rz)
    bestks[[nlt]][[as.character(oka)]] = oteste
  }
}
 
 
#faz uma figura com os resultados
bestks.tb = as.data.frame(bestks.tb)
rownames(bestks.tb) = paste(bestks.tb$FUN,bestks.tb$K)
bestks.tb = bestks.tb[order(as.numeric(bestks.tb$K)),]
#bestks.tb = na.omit(bestks.tb)
fn = paste(pathparafiguras,"/otimoNgruposKmeansTest.pdf",sep="")
vl = is.na(bestks.tb$K)
sum(vl)
bestks.tb = bestks.tb[!vl,]
pdf(fn,width=8.5,height=10)
par(mfrow=c(3,2),mar=c(5,4,3,2))
for(acol in 4:6) {
  clnam = colnames(bestks.tb)[acol]
plot(bestks.tb[!is.na(bestks.tb$FUN),c(2,acol)],lwd=1,type='n',main='Acertos significativos',ylab=paste('Taxa de acerto',clnam))
ff = unique(bestks.tb$FUN)
ff = ff[!is.na(ff)]
ltys = c("solid","dotted","dashed","dotdashed")
pchscol = c("red","blue","green","purple") 
for(f in 1:length(ff)) {
  vll = bestks.tb$FUN==ff[f] & !is.na(bestks.tb[,acol])
  points(bestks.tb[vll,c(2,acol)],lwd=1,col=pchscol[f],type='l')
  points(bestks.tb[vll,c(2,acol)],pch=21,bg=pchscol[f])
}
legend("bottomleft",legend=ff,pch=21,pt.bg=pchscol,cex=0.8,bg='white')
abline(h=0.95,lty='dotted')
plot(bestks.tb[!is.na(bestks.tb$FUN),c(2,acol+3)],lwd=1,type='n',main='Todos os acertos',ylab=paste('Taxa de acerto',clnam))
for(f in 1:length(ff)) {
  vll = bestks.tb$FUN==ff[f] & !is.na(bestks.tb[,acol+3])
  points(bestks.tb[vll,c(2,acol+3)],lwd=1,col=pchscol[f],type='l')
  points(bestks.tb[vll,c(2,acol+3)],pch=21,bg=pchscol[f])
}
legend("bottomleft",legend=ff,pch=21,pt.bg=pchscol,cex=0.8,bg='white')
abline(h=0.95,lty='dotted')
}
dev.off()

############################
#A figura acima permite definir grupos morfológicos segundo algum 
#critério de análise de sensibilidade. O objetivo não é definir 
#espécies, apenas definir subpopulações de dados ('latent populations' em Hennig 2013),
#grupos morfológicos com algum suporte nos dados. O significado biológico desses grupos precisa ser interpretado à luz de outras evidências.
#Melhor portanto pensar: qual o k me dá segurança para falar alguma coisa e interpretá-lo? Isso deve ser feito num bate e volta, entre definir grupos, entender como se diferenciam, como eles se relacionam aos seus pre-conceitos (a priori seus, a priori pelo nome das etiquetas), qual o nível de simpatria. Etc. Com isso na mão você deve tomar decisões sobre quantas espécies e como as amostras devem ser classificadas.
#Exemplo da tomada de decisões com base na figura gerada (análise de sensibilidade):
#a) eu quero K grupos que tem alta taxa de predição em modelos de validação cruzada (loo) por modelos (SVM, Bayes e/ou LDA)
#a.1) no exemplo, quando K=8 e o modelo de gerar grupos é Kmeans a taxa de acerto total é > 95% para SVM, BAYES e LDA. K maiores do que isso tem taxa >95% apenas em Kmeans+LDA (k=11 ou 12). 
#a.2) K=4 é outro agrupamento que tem relativo alto supporte em vários métodos (LDA,SVM,BAYES) e agrupamentos distintos (specc,mclust,kmeans)
 
#eu posso portanto explorar essas três hipóteses (k=8, k=12, k=4)
#Se eu quero uma abordagem splitter ou se minhas idéias iniciais são que meu grupo tem em torno de 11 ou 12 espécies, importante gerar e interpretar o que significa k=12. 
#os diferentes ks explorados aqui permitirao visualizar a hierarquia das similaridades entre os grupos
 
#primeiro vamos obter as classificações para cada K
#depois vamos plotar essas classificações em espaco discriminante
#a figura acima sugere que com K means, nosso poder preditivo é maior para todos esses Ks, então vamos ficar com esse modelo por simplicidade (você pode, se quiser, comparar as classificacoes de diferentes modelos e usar isso para ajudar na classificacao de incertezas)
 
#reclassifica e identifica incertezas para esses Ks
quaisks=c(4,8,12)
 
#cria um objeto para salvar resultados
resd = dado.grupos[,1:2]
#refaz o agrupamento para esse valor de K (número de grupos)
for(q in 1:length(quaisks)) {
  #com Kmeans
  ttk = kmeans(dado.grupos, center=quaisks[q],iter.max=10000,nstart=1000)
  clk = ttk$cluster
  #com Mclust
  #ttmc = Mclust(dado.grupos,G=quaisks[q])
  #clmc = ttmc$classification
  #tbb = table(KMEANS=clk,MCLUST=clmc)
  #essas classificacoes sao parecidas?
 
  #testa essa classificacao
  dd=dado.teste
  grps = clk
  dd$GROUPING=paste("GRUPO",sprintf("%02s",grps))
  #com lda e loo
  oteste = testa.classificacao(dad=dd,"GROUPING",analises=c("svm","bayes","lda")[3],method=c("loo","perc")[1])
 
  #pega o resultado
  br = oteste$BRUTOS
  #classifica as amostras
  br$CLASSIFICACAO = NA
  #quais amostras tem predicao correta nos grupos
  vl = br$MELHOR.PREDICAO.PROB>=0.95 & br$ACERTO=='acerto' & !is.na(br$ACERTO)
  sum(vl)
  br$CLASSIFICACAO[vl] = "core"
 
  #quais sao acertos nao significativos
  vl = br$MELHOR.PREDICAO.PROB<=0.95 & br$ACERTO=='acerto' & !is.na(br$ACERTO)
  sum(vl)
  br$CLASSIFICACAO[vl] = "incerteza"
  #quais amostras são erros 
  vl = br$ACERTO=='ERRO' & !is.na(br$ACERTO)
  sum(vl)
  #desclassifica os erros
  br$MELHOR.PREDICAO[vl] = NA
  br$MELHOR.PREDICAO.PROB[vl] = NA
  br$CLASSIFICACAO[vl] = "sem classe"
  cls = c("MELHOR.PREDICAO","MELHOR.PREDICAO.PROB","CLASSIFICACAO")
  brr = br[,cls]
  colnames(brr) = paste("K=",quaisks[q],"_",cls,sep="")
  brr = brr[-1,]
  rnn = rownames(brr)
  resd = cbind(resd[rnn,],brr[rnn,])
}
resd = resd[,-c(1:2)]
head(resd)
#salva essas classificacoes
fn = paste(pathparadados,"/classificacaoAposterioriK.csv",sep="")
write.table(resd,file=fn,sep="\t",na="",quote=T,row.names=T)

Compara grupos & Refina classificação final