#gera matrizes de confusao para cada par
#comparando classificacoes a posteriori. Tanto grupos core como grupos com incertezas
fn = "classificacaoaposteriorik.csv"
resd = read.table(file=fn,sep="\t",na.strings=c("NA",""),row.names=1,header=T,as.is=T)
colunasgrupos = seq(1,ncol(resd)-2,by=3)
comparacoes = NULL
tbbs = NULL
for(cl in colunasgrupos) {
resto = colunasgrupos[!colunasgrupos==cl]
for(cl2 in resto) {
comparacao =paste(sort(colnames(resd)[c(cl,cl2)]),collapse="-")
if (!comparacao%in%comparacoes) {
comparacoes = c(comparacoes,comparacao)
tb = table(resd[,cl],resd[,cl2])
tbbs[[comparacao]] = tb
#apenas core
vl = resd[,cl+2]=='core' & !is.na(resd[,cl+2])
vl2 = resd[,cl2+2]=='core' & !is.na(resd[,cl2+2])
vll = vl & vl2
tb = table(resd[vll,cl],resd[vll,cl2])
tbbs[[paste(comparacao,"CORE")]] = tb
}
}
}
fn = paste(pathparafiguras,"/matrizConfusaoEntreKsAposteriori.pdf",sep="")
pdf(file=fn,width=8.5,height=6)
for(t in 1:length(tbbs)) {
print(t)
par(mar=c(2,7,7,2),mfrow=c(1,1),cex.axis=0.8)
tb = tbbs[[t]]
comparacao = names(tbbs)[t]
#define cores do backgroun e texto
bg.cols= tb
txt.cols = tb
#se for 0 é branco
bg.cols[tb==0] = gray(level=1)
txt.cols[tb==0] = gray(level=1)
#se for 1 é cinza claro
bg.cols[tb==1] = gray(level=0.8)
txt.cols[tb==1] = gray(level=0)
#se for entre 1 e 5 é um pouco mais escuro
bg.cols[tb>1 & tb<=5] = gray(level=0.6)
txt.cols[tb>1 & tb<=5] = gray(level=0)
bg.cols[tb>5 & tb<=10] = gray(level=0.4)
txt.cols[tb>5 & tb<=10] = gray(level=1)
#se for > 10 o fundo é preto e o texto branco
bg.cols[tb>10] = gray(level=0)
txt.cols[tb>10] = gray(level=1)
#one for 0 coloca NA, isso não é plotado
tb[tb==0] = NA
plotamatriz(tb,bg.cols=bg.cols,txt.cols=txt.cols,valcex=0.8, cexaxis = 0.8)
mtext(side=3,text=comparacao,adj=0,line=6,cex=0.8)
}
dev.off()
rm(list=ls())
source("funcoesNecessarias.R")
#para salvar resultados
pathparadados = "dadosParaAnalises"
pathparafiguras = "figuras"
pathparatabelas = "tabelas"
#le dados
arq = "dadosParaAnalises/dadosMORFOnmdsvarsel.csv"
dados.morfo = read.table(file=arq,sep='\t',as.is=F,row.names = 1,header=T)
#le os dados e gera matrizes de contigencia (confusao)
fn = paste(pathparadados,"/classificacaoAposterioriK.csv",sep="")
resd = read.table(file=fn,sep="\t",na.strings=c("NA",""),row.names=1,header=T,as.is=T)
#reclassificando:
#1. Encontrar congruencia entre classificacoes diferentes (K diferentes) pode indicar amostras que sempre vao juntas e amostras que "pulam" de grupo. Dessa forma seria possível:
#a) definir como modelo amostras congruentemente classificadas (CORE);
#b) reclassificar as demais amostras com base nesse modelo (refinar ate estabilizar)
#c) checar novamente a congruencia e classificar as amostras em grupo
#c1)"core" aqueles subgrupos que sempre vao juntos
#c2)"incertezas" amostras ou grupos de amostras que nem sempre vao juntos
#vamos fazer isso com duas classificacoes feitas
#a posteriori (k=4 e k=8)
cln1 = "K.4_MELHOR.PREDICAO"
cln2 = "K.8_MELHOR.PREDICAO"
cln1p = "K.4_MELHOR.PREDICAO.PROB"
cln2p = "K.8_MELHOR.PREDICAO.PROB"
vl1 = resd[,cln1p]>0.95 & !is.na(resd[,cln1p])
vl2 = resd[,cln2p]>0.95 & !is.na(resd[,cln2p])
#quais tem alta predicao em ambas classificacoes
vl = vl1 & vl2
#qual a matriz de confusao para essas amostras
#grupos sao vetores com nomes de linha
grupo1 = resd[vl,cln1]
names(grupo1) = rownames(resd)[vl]
grupo2 = resd[vl,cln2]
names(grupo2) = rownames(resd)[vl]
tb = table(grupo1,grupo2)
plotamatriz(tb)
#neste caso são dois grupos hierárquicos de uma mesma metodologia (kmeans).
#note que a hierarquia não é perfeita. Membros do grupo 1 e 6 de k=8 tem amostras em quase todos os grupos de k=4. Mas grupo 1 e 6 tem uma maioria das amostras concentrada num grupo de k=4.
#para este exemplo, portanto, vamos reclassificar as amostras de k=8, deixando apenas aquelas que aparecem juntas em k=4, sem criar novos grupos. As demais amostras serão classificadas como incertezas após este refinamento
#refina isso gerando modelos com os grupos core
#e predizendo as demais
#veja o codigo da funcao para entender o que isso esta fazendo
grupo1 = resd[,cln1]
grupo1.prob = resd[,cln1p]
grupo2 = resd[,cln2]
grupo2.prob = resd[,cln2p]
rn = rownames(resd)
names(grupo2) = rn
names(grupo1) = rn
names(grupo1.prob) = rn
names(grupo2.prob) = rn
res = refinagrupos(grupo1,grupo1.prob,grupo2,grupo2.prob,dados.morfo)
#veja o que isso gerou
#tabulando novamente apenas para dados core
#o match deve ser perfeito
tbb = table(res$GRUPO1.ORG,res$NOVACLASSE)
plotamatriz(tbb)
#plota apenas core (neste caso o match deve ser total)
vl = res$CATEGORIA=='core'
tbb = table(res$GRUPO1.ORG[vl],res$NOVACLASSE[vl])
plotamatriz(tbb)
#note que esses grupos core (k=8) batem com apenas 1 grupo de k=4, ou seja essas amostras vao juntas nao importa a escala
#faz um gráfico morfológico (espaço LDA)
#mostrando grupos e incertezas
#pega eixos discriminantes
gp = res$NOVACLASSE
mld = lda(gp ~ ., dados.morfo)
#pega os eixos discriminantes
pd = predict(mld)
tt = pd$x[,1:3]
#define cores para k grupos
cores = rainbow(n=8)
#define simbolos para k gps
pchs.core = c(21,22,23,24,25,21,22,23)
pchs.inc = c(1,0,5,2,6,1,0,5)
#define tamanhos segundo as probabilidades
pr = round(res$PROBABILIDADE,2)
#mantem tamanho 0.5 para a maioria (que tem predicao alta)
vl = pr>0.95
cexx = pr
cexx[vl] = 0.8
#as demais coloca tanto maior quanto menor a probabilidade
cexx[!vl] = 0.6+(1-cexx[!vl])*3
#plota o espaco lda e 3 dimensoes (duas figuras)
cexx = cexx+(cexx*0.2)
fn = paste(pathparafiguras,"/ldaGrupoFinal.pdf",sep="")
pdf(fn,width=8.5,height=10)
par(mar=c(5,5,2,1),mgp=c(2.5,0.5,0),tck=-0.015,cex.axis=0.8)
mt = matrix(c(1,1,1,2,3,3,3,4),nrow=2,ncol=4,byrow = T)
lm =layout(mt)
LDA1=tt[,1]
LDA2=tt[,2]
#plota vazio
plot(LDA1,LDA2,type='n')
gp = as.factor(gp)
vl = res$CATEGORIA=='core'
points(LDA1[vl],LDA2[vl],pch=pchs.core[as.numeric(as.factor(res$NOVACLASSE[vl]))],bg=cores[as.numeric(as.factor(res$NOVACLASSE[vl]))],cex=cexx[vl])
#incertezas
vl = !vl
points(LDA1[vl],LDA2[vl],pch=pchs.inc[as.numeric(as.factor(res$NOVACLASSE[vl]))],col=cores[as.numeric(as.factor(res$NOVACLASSE[vl]))],cex=cexx[vl])
#plota grupo menor
lvs = levels(res$GRUPO1.ORG)
for(lv in lvs) {
vl = res$GRUPO1.ORG==lv & !is.na(res$GRUPO1.ORG) & res$CATEGORIA=='core'
dld = tt[vl,1:2]
xy = chull(dld)
xy <- c(xy, xy[1])
lines(dld[xy,],lty='dotted')
}
plot.new()
par(mar=c(3,0,3,1))
legend("left",legend=c(levels(gp),paste(levels(gp),'incertezas')),pch=c(pchs.core,pchs.inc),col=c(rep('black',nlevels(gp)),cores),pt.bg=cores,pt.cex=1,ncol=1,cex=1,y.intersp = 1.2,bty='n')
legend("bottomleft",legend=c("grupos k=4"),lwd=2,lty='dotted',bty='n',cex=1)
LDA1=tt[,1]
LDA3=tt[,3]
#plota vazio
par(mar=c(5,5,2,1))
plot(LDA1,LDA3,type='n')
vl = res$CATEGORIA=='core'
points(LDA1[vl],LDA3[vl],pch=pchs.core[as.numeric(as.factor(res$NOVACLASSE[vl]))],bg=cores[as.numeric(as.factor(res$NOVACLASSE[vl]))],cex=cexx[vl])
#incertezas
vl = !vl
points(LDA1[vl],LDA3[vl],pch=pchs.inc[as.numeric(as.factor(res$NOVACLASSE[vl]))],col=cores[as.numeric(as.factor(res$NOVACLASSE[vl]))],cex=cexx[vl])
gp = as.factor(gp)
lvs = levels(res$GRUPO1.ORG)
for(lv in lvs) {
vl = res$GRUPO1.ORG==lv & !is.na(res$GRUPO1.ORG) & res$CATEGORIA=='core'
dld = tt[vl,c(1,3)]
xy = chull(dld)
xy <- c(xy, xy[1])
lines(dld[xy,],lty='dotted')
}
plot.new()
par(mar=c(3,0,3,1))
legend("left",legend=c(levels(gp),paste(levels(gp),'incertezas')),pch=c(pchs.core,pchs.inc),col=c(rep('black',nlevels(gp)),cores),pt.bg=cores,pt.cex=1,ncol=1,cex=1,y.intersp = 1.2,bty='n')
legend("bottomleft",legend=c("grupos k=4"),lwd=2,lty='dotted',bty='n',cex=1)
par(mar=c(2,7,7,2))
mt = matrix(c(1,2),nrow=2,ncol=1,byrow = T)
layout(mt)
tbb = table(res$GRUPO1.ORG,res$NOVACLASSE)
plotamatriz(tbb)
mtext(side=3,text="todas",line=3,adj=-0.1,font=2,cex=1.5)
vl = res$CATEGORIA=='core'
tbb = table(res$GRUPO1.ORG[vl],res$NOVACLASSE[vl])
plotamatriz(tbb)
mtext(side=3,text="core",line=3,adj=-0.1,font=2,cex=1.5)
dev.off()
#salva essa classificacao final
fn = paste(pathparadados,"/classificacaoAposterioriFINAL.csv",sep="")
write.table(res,file=fn,sep="\t",na="",quote=T,row.names=T)