Escolhendo um valor de K para o Kmeans usando o critério de Calinski

No último post, aqui, a gente viu que quando agrupamos espécies em conjuntos, que tenham algum significado para a gente, existe coisas interessantes que podem ser feitas, como achar espécies que identifiquem aquele conjunto.

Mas isso é um tipo de classificação supervisionada, já que nos classificamos os locais de coleta, a priori, cada ponto está num grupo que nós definimos. Agora suponha que achemos que existem grupos para um conjunto de amostras, mas não sabemos a que grupo cada amostra está associada? E agora?

Existe um algorítimo chamado kmeans, que podemos usar. Mas o foco agora não é explicar como ele funciona, mas como escolher quantos grupos existem nos dados. Vamos tentar pensar num exemplo para ser mais claro. Pense que temos duas espécies só (ou se quiser, pense em dois atributos de uma flor, tamanho de pétalas e sépalas, similar aos dados do iris usado pelo Fisher), x e y, e podemos representá-las em um gráfico 2d.

01

Então temos quatro grupos, que podemos identificar visualmente na figura certo? Cada um dos grupos vem de distribuições normal com média e desvios diferentes, tanto para x quanto para y.

Note também que o preto e o verde estão um pouquinho mais embolados entre si (nos que quisemos assim) em relação ao vermelho e azul. Mas aqui essencialmente temos 120 amostras, e cada uma pertence a um dos 4 grupos que definimos no nosso exemplo, ao simular os dados.

Até aqui tudo bem, mas imagine, que não soubéssemos a cor de cada ponto, não soubéssemos a que grupo cada um pertence, e agora? Veja como é isso abaixo.

02

Aqui é um exemplo trivial, então só no olhômetro da para ver que existem grupos nessas 120 amostras, alias quatro grupos. Agora podemos usar o k-means para atribuir um grupo a cada uma dessas.

Então nos mandamos o kmeans classificar os dados em 4 grupos.

1
ajuste<-kmeans(exemplo,4)

Simples assim, e podemos comparar o resultado com os grupos originais.

03

Veja que a classificação erra um pouco, principalmente nos pontos verdes que estão misturados no pretos (lembrando que uma cor é um grupo, a gente troca as cores, porque o que interessa é pertencer ao grupo certo, não ter a cor certa), mas de modo geral, o kmeans recupera bem o grupo original de cada amostra. Claro que ele tem uma série de pressupostos a serem seguidos, mas isso fica para outro post.

Agora problema é que eu coloquei para separar em 4 grupos, porque eu quis assim. E isso é um problema, como a gente vai saber quantos grupos um conjunto de dados tem? Veja que assim como eu usei 4, eu poderia ter usado outras quantidades de grupos.

04

Lembrado, que quando realmente vamos usar o kmeans, a gente não sabe quantos grupos é a melhor opção. Ai entra o critério de Calinski, é uma medida de o quão bom é o ajuste de um determinado número de grupos, de forma que podemos descobrir quantos grupos usar. Podemos discernir qual dos quadros acima é mais adequado, usando o critério e todos e vendo quando ele é maximizado.

Ele está implementado no pacote vegan do R, na função cascadeKM, e é usado assim:

1
ajuste<-cascadeKM(exemplo,2,10,criterion = 'calinski')

A gente fala a nossa tabela de dados, o número mínimo de grupos, o número máximo e qual critério usar. Dai para cada valor de k(número de grupos), o cascadeKM vai rodar o kmeans, calcular o valor de calinski e dizer a qual grupo cada amostra esta nessa classificação.

Basicamente ele calcula o seguinte:

calinski=\frac{SSB/(k-1)}{SSW/(n-k)}

onde n é o número de amostras e k é o número de classes que estamos selecionando (que vai ser 2,3,4,…10, veja que colocamos 2 e 10 como argumentos). SSW é a soma dos quadrados dentro dos clusters e SSB é a soma dos quadrados entre clusters, sendo análogo a anova (que lembrem que F, quanto mais, mais quer dizer que existem grupo, veja aqui), ou seja, a gente vai procurar o maior valor de F.

Por padrão, o cascadeKM faz o seguinte plot

05

Que é a qual grupo cada ponto foi atribuído e uma figura com os valores de calinski e o número de k. Mas ele está invertido, então podemos também extrair os valores e produzir uma figura com o valor de k no eixo x e o resultado a métrica de calinski no eixo y.

06

Veja que o maior valor está exatamente no 4, que é como geramos os dados. Porém nem tudo é tão simples, esse critério vai se comportar melhor se os dados estiverem em uma região especifica do n-espaço dos dados, de preferencia com variação homogenia, além de outros casos. Mas de qualquer forma, ele pode ser informativo, sendo que podemos avaliar o resultado graficamente para ver se faz sentido, e ele é bastante simples de entender e calcular.

Bem é isso ai, como estava olhando essas coisas, resolvi aproveitar e deixar tudo que estava testando como post, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter R. Minchin, R. B. O’Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens and Helene Wagner (2016). vegan: Community Ecology Package. R package version 2.3-5. https://CRAN.R-project.org/package=vegan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
##Pacotes utilizados
library(vegan)
 
##Gerando um exemplo
x<-c(1,1,2,2)
y<-c(1,2,1,2)
 
set.seed(123)
exemplo<-data.frame(
    x=c(rnorm(30,x[1],0.3),rnorm(30,x[2],0.2),rnorm(30,x[3],0.3),rnorm(30,x[4],0.2)),
    y=c(rnorm(30,y[1],0.2),rnorm(30,y[2],0.2),rnorm(30,y[3],0.2),rnorm(30,y[4],0.2))
)
grupo<-factor(rep(letters[1:4],each=30))
 
##Exemplo gerado
plot(exemplo$x,exemplo$y,col=as.numeric(grupo),pch=19,frame=F,xlab="eixo x",ylab="eixo y")
legend("topright",legend=c(levels(grupo)),col=unique(as.numeric(grupo)),pch=19,bty="n")
 
##Dados não classificados
plot(exemplo$x,exemplo$y,pch=19,frame=F,xlab="eixo x",ylab="eixo y")
 
ajuste<-kmeans(exemplo,4)
 
##Resultado vc original
par(mfrow=c(2,1))
plot(exemplo$x,exemplo$y,col=as.numeric(grupo),pch=19,xlab="eixo x",ylab="eixo y")
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y")
 
##Comparando varios valores de k numa única figura
par(mfrow=c(2,2),mar=c(4,4,2,2))
ajuste<-kmeans(exemplo,2)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=2")
ajuste<-kmeans(exemplo,3)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=3")
ajuste<-kmeans(exemplo,6)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=6")
ajuste<-kmeans(exemplo,8)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=8")
 
##
ajuste<-cascadeKM(exemplo,2,10,criterion = 'calinski')
 
##Plot original
plot(ajuste)
 
##Fazendo nosso próprio plot
plot(2:10,ajuste$results[2,],type="b",pch=19)
abline(v=which.max(ajuste$results[2,])+1,lty=3,lwd=3,col="red")

Espécies indicadores de Agrupamento (indicator Value ou indval)

Uma atividade comum em ecologia, assim como em muitas áreas do conhecimento, é agrupar coisas.
Tipos de ambientes, grupos de espécies, seja la o que for, baseado em atributos, que podem ser características morfológicas ou abundancia de espécies.

Vamos pensar num exemplo, imagine que coletamos 30 amostras, fomos em 30 locais, 15 dos quais são reservas para conservação e outras 15 que são áreas com atividades extrativistas.

Então temos um conjuntos de dados mais ou menos assim:

> comunidade Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 Local 1 0 2 0 0 4 2 Local 2 2 1 0 0 4 1 Local 3 1 1 0 0 2 1 Local 4 2 1 0 0 2 0 Local 5 3 0 0 0 1 1 Local 6 0 0 0 0 4 4 Local 7 1 1 0 0 1 0 Local 8 2 0 0 0 0 0 Local 9 1 0 0 1 5 0 Local 10 1 0 0 0 3 2 Local 11 3 0 0 0 1 2 Local 12 1 0 0 2 2 3 Local 13 1 0 0 1 5 2 Local 14 1 0 0 1 2 2 Local 15 0 0 0 0 1 1 Local 16 2 0 2 0 0 0 Local 17 0 0 0 2 0 1 Local 18 0 0 1 1 0 1 Local 19 0 0 1 2 0 2 Local 20 3 1 3 1 0 0 Local 21 2 0 2 0 0 0 Local 22 1 0 3 2 0 0 Local 23 1 1 3 0 0 0 Local 24 4 0 3 1 0 0 Local 25 1 0 2 1 0 1 Local 26 1 0 3 1 0 0 Local 27 1 0 2 1 0 0 Local 28 1 1 3 1 0 0 Local 29 0 1 0 4 0 0 Local 30 0 0 2 1 0 1

Primeiro, para vermos se nossas classes são relevantes, se realmente existe uma diferença, podemos primeiro tentar visualizar esses dados. Uma maneria comum é usar alguma métrica de similaridade, ou dissimilaridade, que da, quase sempre, na mesma. Então a gente calcula alguma distancia, como a de Bray-Curtis para os locais e podemos fazer uma árvore baseado nas distâncias.

01

Agora veja que cada local pertence a uma classe, um local de origem no nosso caso, então por exemplo, o local 1 e 2 são da reserva conservada, enquanto locais 15 a 30 são da parte extrativista. Então podemos substituir esses nomes na árvore.

02

Veja que parece que temos um grupo de perturbados e outro grupo de conservados se formando, mas temos alguns locais que não são tão óbvios, pequenos grupos meio mistos, meio bagunçados.

Uma outra forma de visualização é com algum tipo de redução dimensional, sendo os mais conhecidos o escalonamento multidimensional e análises baseadas em vetores eigen. Vamos realizar uma análise de componentes principais aqui, para ver como fica essa comunidade.

03

Aqui, geramos vários eixos, que representam a similaridade entre os locais, sendo que quanto mais próximo dois pontos estão, mais similares eles são. Além disso, já trocamos o nome do local pelo grupo que ele pertence. Existem várias coisas que deveríamos observar, antes de olhar dessa figura, mas vamos pular essa parte. Veja que temos algo como dois grupos até bem definidos aqui, se passarmos uma linha vertical no primeiro eixo, vemos fácil essa separação, que é fácil de ver porque simulados os dados dessa forma, lembre-se disso. Mas veja que realmente, visualmente, tanto usando o dendograma da árvore como a análise de componentes principais vemos que realmente esse grupos parecem existir.

Mas agora vem a pergunta, beleza, existem grupos, mas se pegarmos uma nova amostra, sem saber o grupo dela, como podemos saber qual o grupo dela, sem saber as coordenadas de coleta?
Uma primeira observação a fazer é se alguma espécie é comum em uma área e incomum na outra. E isso pode ser feito calculando o valor de indicação proposto por Dufrene-Legendre.

Existem alguns pacotes no R que tem esse valor já implementado em funções fáceis de usar, entre eles existem o labdsv e o indicspecies.

No labdsv, isso é bem fácil.

Precisamos dos dados de comunidade, em uma tabela, como no inicio do post e um vetor com o agrupamento, a classe de cada local, como por exemplo o fator abaixo.

> agrupamento [1] Conservado Conservado Conservado Conservado Conservado Conservado Conservado [8] Conservado Conservado Conservado Conservado Conservado Conservado Conservado [15] Conservado Pertubado Pertubado Pertubado Pertubado Pertubado Pertubado [22] Pertubado Pertubado Pertubado Pertubado Pertubado Pertubado Pertubado [29] Pertubado Pertubado Levels: Conservado Pertubado

Com essas duas informações, comunidade e classes, podemos calcular o “Indicator Value” e realizar um teste de hipóteses para avaliar sua significância, se realmente a espécie pode ser uma boa opção para diferencias as classes, baseado na sua ocorrência e abundância.

> indval(comunidade,agrupamento) $relfrq Conservado Pertubado Sp 1 0.8000000 0.6666667 Sp 2 0.3333333 0.2666667 Sp 3 0.0000000 0.8666667 Sp 4 0.2666667 0.8000000 Sp 5 0.9333333 0.0000000 Sp 6 0.7333333 0.3333333 $relabu Conservado Pertubado Sp 1 0.5277778 0.4722222 Sp 2 0.6000000 0.4000000 Sp 3 0.0000000 1.0000000 Sp 4 0.2173913 0.7826087 Sp 5 1.0000000 0.0000000 Sp 6 0.7777778 0.2222222 $indval Conservado Pertubado Sp 1 0.42222222 0.31481481 Sp 2 0.20000000 0.10666667 Sp 3 0.00000000 0.86666667 Sp 4 0.05797101 0.62608696 Sp 5 0.93333333 0.00000000 Sp 6 0.57037037 0.07407407 $maxcls Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 1 1 2 2 1 1 $indcls Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 0.4222222 0.2000000 0.8666667 0.6260870 0.9333333 0.5703704 $pval Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 0.618 0.732 0.001 0.008 0.001 0.006 $error [1] 0 attr(,"class") [1] "indval"

A saída da função é meio longa, mas duas partes nos importam mais, primeiro o $indval.
O indival é o indicator value de cada espécie, para cada nivel da categoria, veja que temos uma categoria chamada agrupamento, e essa categoria tem dois níveis, conservado e pertubado.

Então, para cada espécie, para cada nível, a gente calcula um valor, que diz o quanto aquela espécie é relevante como classificadora daquele nível. Ou seja, se eu ver essa espécie aqui implica que esse local é de tal categoria. Mas calma la, vamos ver essa tabela melhor.

$indval Conservado Pertubado Sp 1 0.42222222 0.31481481

Para a espécie 1, ela tem um valor tanto para a área conservado, como para a área perturbada. Ou seja, apesar dela ter um sinal indicador, ela não me resolve nada, já que ela indica todos os níveis.

Agora olha que legal, a espécie 3.

$indval Conservado Pertubado Sp 1 0.42222222 0.31481481 Sp 2 0.20000000 0.10666667 Sp 3 0.00000000 0.86666667

Ela tem uma valor alto no perturbado, e além disso tem zero na área Conservado, então ela deve estar ocorrendo somente na área perturbada, e se a gente voltar e olhar na tabela de dados, veremos que é isso mesmo, claro que aqui, como temos 6 especies, 30 locais, isso é fácil de ver, mas isso começa a complicar, basta aumentar a quantidade de locais e espécies, que olhar a tabela se torna uma atividade difícil, e essa métrica se torna muito útil.

Agora veja que temos um resumo desses valores, que seria útil no caso de termos muitas classes

$indcls Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 0.4222222 0.2000000 0.8666667 0.6260870 0.9333333 0.5703704

O $indcls é então somente o maior valor de cada linha naquela matriz $indval, e temos também um teste de hipoteses, para saber quem é releventa para fazer uma classificação, poe sabemos que a espécie 1 ajudaria pouco a classificar o local, a indicar o nível que pertence o local, mas o espécie 3 ajuda muito.

$pval Sp 1 Sp 2 Sp 3 Sp 4 Sp 5 Sp 6 0.618 0.732 0.001 0.008 0.001 0.006

E vemos que a espécie 1 tem não tem uma significância, se fizermos um corte no 0.05, mas a espécie 3 sim, tem uma significância.

Seguindo a mesma ideia, mas com algumas modificações na ideia do indval, temos a análise feita pelo pacote indicspecies, que na prática tem uma saída com menos coisas para interpretar, mas ainda sim similar.

> indval <- multipatt(comunidade, agrupamento,control = how(nperm=5000)) > summary(indval) Multilevel pattern analysis --------------------------- Association function: IndVal.g Significance level (alpha): 0.05 Total number of species: 6 Selected number of species: 4 Number of species associated to 1 group: 4 List of species associated to each combination: Group Conservado #sps. 2 stat p.value Sp 5 0.966 0.0002 *** Sp 6 0.755 0.0146 * Group Pertubado #sps. 2 stat p.value Sp 3 0.931 0.0002 *** Sp 4 0.791 0.0046 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Temos por nível do agrupamento, quais espécies são significativamente relevantes, são omitidas as espécies irrelevantes, para a classificação do local. Veja que o stat, é diferente dos valores de indval calculados pela função do pacote labdsv.

Bem é isso ai, uma forma legal de identificar quem pode ser util para classificação de um ambiente, mas exite também a possibilidade de usar árvores de classificação para essa atividade. O script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Dufrene, M. and Legendre, P. 1997. Species assemblages and indicator species: the need for a flexible asymmetrical approach. Ecol. Monogr. 67(3):345-366.

De Cáceres, M. and P. Legendre 2009. Associations between species and groups of sites: indices and statistical inference. Ecology 90: 3566-3574.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
library(labdsv)
library(indicspecies)
 
set.seed(123)
agrupamento <- factor(rep(c("Conservado","Pertubado"),each=15))
 
n_amostras<-30
n_sp<-6
 
sp1<-rpois(n_amostras,1)
sp2<-rpois(n_amostras,0.5)
sp3<-rpois(n_amostras,(as.numeric(agrupamento)-1)*2)
sp4<-rpois(n_amostras,(as.numeric(agrupamento)-0.65))
sp5<-rpois(n_amostras,(as.numeric(agrupamento)-2)*-2)
sp6<-rpois(n_amostras,rev(as.numeric(agrupamento)-0.65))
 
comunidade<-as.data.frame(matrix(c(sp1,sp2,sp3,sp4,sp5,sp6),ncol=n_sp,nrow=n_amostras,byrow=F,
                                 dimnames=list(paste("Local",1:n_amostras),paste("Sp",1:n_sp))))
comunidade
dis.comunidade <- dsvdis(comunidade,'bray/curtis')
cluster<-hclust(dis.comunidade,method="complete")
 
##
jpeg("01.jpg")
plot(cluster)
dev.off()
##
jpeg("02.jpg")
plot(cluster,labels=agrupamento)
dev.off()
 
pca_comunidade<-princomp(comunidade)
 
jpeg("03.jpg")
plot(pca_comunidade$scores[,1],pca_comunidade$scores[,2],type="n",xlab="Eixo 1",ylab="Eixo2",frame=F)
text(pca_comunidade$scores[,1],spca_comunidade$scores[,2],agrupamento)
abline(v=0,col="red",lty=3,lwd=2)
dev.off()
 
jpeg("04.jpg")
biplot(pca_comunidade)
dev.off()
 
agrupamento
indval(comunidade,agrupamento)
 
 
indval <- multipatt(comunidade, agrupamento,control = how(nperm=5000))
summary(indval)