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)

Determinando a distribuição espacial a partir da Distribuição do vizinho mais próximo

Quando a gente estava medindo as distâncias ao vizinho mais próximo aqui, se a gente olhar a distribuição desses, a gente vai ver algo interessante aparecer.

Primeiro, vamos transformar em funções aquelas operações.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
#Calcula distancias do vizinho mais próximo
calcula_r<- function(x,y){
    r <- vector()
    nn <- vector()
    d <- vector()
    for (i in 1:100) {
        d <- distancia(x[i],y[i],x,y)
        r[i] <- min(d[-i])
    }
    return(r)
}

O que nos interessa, é calcular a média das distâncias ao vizinho mais próximo.

1
2
3
4
5
6
7
n<-1000
dist_vmp<-vector()
 
for(i in 1:n){
    r<-calcula_r(runif(100),runif(100))
    dist_vmp[i]<-mean(r)
}

Fazemos isso muitas vezes, uma mil, e podemos dar uma olhada em como fica a distribuição dessas mil médias. Veja que estamos gerando uma distribuição de pontos ao acaso no espaço, e para cada um desses conjuntos, a gente calcula a média dos vizinhos mais próximos.

01

Veja que é assim que fica a média de vizinhos mais próximos, quando os pontos estão distribuídos ao acaso no espaço. Veja que se a gente pegar outro exemplo, e vamos ilustrar os pontos aqui.

02

Em relação a nossa distribuição, veja onde fica o média de vizinhos mais próximos.

03

Bem, até aqui, nada de novo, a gente sempre está gerando os pontos ao acaso, então temos essa distribuição média. Agora vamos considerar um segundo caso.

04

Agora nos temos os pontos agregados em alguns pontos. Agregações na natureza ocorrem por muitos motivos, pode ser mais seguro ficar perto de outros pontinhos, ali pode ser o local onde os melhores recursos estão estocados, enfim várias coisas ponde acontecer. Mas isso tem algum impacto na média do vizinhos mais próximo? Bem, todo mundo vai tender a estar mais pertinho, e assim teremos uma média mais baixa. Para esse exemplo, podemos calcular a média e comparar com a distribuição dos mil casos aleatórios.

05

E o que acontece é exatamente isso, o valor de média fica bem menor, bem antes da distribuição, já que todo mundo mais pertinho que o acaso. Porque algum motivo atrai todo os pontinho para o mesmo lugar.

Agora existe outra possibilidade, que é algo como, melhor sozinho do que acompanhado, ou seja, pense num caso de competição extrema, então a gente tenta ficar o mais longe possível de outros competidores.

07

Que gera um padrão mais ou menos como o acima. Agora o que vai acontecer com a média de distância do vizinho mais próximo? Agora ela fica um valor alto, porque sempre temos valores altos.

08

O que estamos vendo aqui? Baseado na média da distancia do vizinho mais próximo, podemo ter uma noção da distribuição espacial dos pontos, sabendo se ela é aleatória ou não, e não sendo aleatória, por esse valor, temos ainda uma noção se tudo está agregado ou espaçado. Sendo que existe um gradiente entre agregado e espaço, e o aleatório está no meio desse gradiente.

09

Teoricamente, a gente nem precisa calcular a distribuição aleatória dos pontos, porque temos uma expectativa para esse valor, que é

E(r)=\frac{\sqrt{p}}{2}

Onde r é a média das distancias do vizinho mais próximo, E é de esperança, a média para distribuições e p é a densidade de pontos, quantos pontos temos por unidade de espaço.

Se a gente olhar a média da nossa simulação

> mean(dist_vmp) [1] 0.05220953

E sabendo que temos 100 pontos, para uma área de 10000 cm2, área aqui é uma unidade meio vaga ja que estamos numa simulação.

> E_r<-sqrt(100/10000)/2 > E_r [1] 0.05

Mas chegamos a um valor bem próximo de média. E como temos uma expectativa, isso facilita a tentar interpretar esse valor como uma razão, se dividirmos o valor que encontramos por esse valor esperado, temos que:

Nos casos de distribuição aleatória, ele estará próximo de 1

> r_ale/E_r [1] 1.111665

Nos casos de agregamento, ele estára abaixo de 1

> r_agre/E_r [1] 0.4468558

E no caso de distribuição uniforme, ele será maior que 1

> r_uni/E_r [1] 1.810889

Mas apesar de tudo essa técnica leva a um outro problema, que é o seguinte. Para conseguirmos ter uma ideia do que é esperado para o local, precisamos saber a densidade de pontos do local, mas isso a gente só sabe com certeza, se tivermos certeza que estamos vendo todos os pontos do local, e como já vimos varias vezes, nossas contagens são normalmente tudo menos precisa. Logo isso é um problema para ter um valor esperado para r, o que dificulta bastante o uso desse método. Mas na teoria ele funciona bem.

Bem é isso ai, brincando mais um pouquinho com as distancias do vizinho mais próximo, 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:
Michael J Crawley 2007 – R Book 2nd. Editora Wiley 1072pp .

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
#Calcula distancias do vizinho mais próximo
calcula_r<- function(x,y){
    r <- vector()
    nn <- vector()
    d <- vector()
    for (i in 1:100) {
        d <- distancia(x[i],y[i],x,y)
        r[i] <- min(d[-i])
    }
    return(r)
}
 
##Simular média do vizinho mais próximo
n<-1000
dist_vmp<-vector()
 
for(i in 1:n){
    r<-calcula_r(runif(100),runif(100))
    dist_vmp[i]<-mean(r)
}
 
##Resultado da simulação
hist(dist_vmp,main="Distribuição de médias do vizinho mais próximo")
 
##Avaliando alguns casos:
##Caso aléatorio
x<-runif(100)
y<-runif(100)
 
plot(x,y)
 
r_ale<-mean(calcula_r(x,y))
hist(dist_vmp)
abline(v=r_ale,lwd=2,col="blue",lty=2)
legend("topright",lwd=2,col="blue",lty=2,legend="Aleatório",bty="n")
 
##Caso Agregado
x<-c(rnorm(50,0.25,0.05),rnorm(50,0.6,0.1))
y<-c(rnorm(50,0.25,0.1),rnorm(50,0.7,0.05))
 
plot(x,y)
 
r_agre<-mean(calcula_r(x,y))
hist(dist_vmp,xlim=c(min(dist_vmp,r_agre),max(dist_vmp)))
abline(v=r_agre,lwd=2,col="red",lty=2)
legend("top",lwd=2,col="red",lty=2,legend="Agregado",bty="n")
 
##Caso Uniforme
x<-rep(seq(0.01,0.99,length=10),10)+runif(100,-0.002,0.05)
y<-rep(seq(0.01,0.99,length=10),each=10)+runif(100,-0.002,0.05)
 
plot(x,y)
 
r_uni<-mean(calcula_r(x,y))
hist(dist_vmp,xlim=c(min(dist_vmp),max(dist_vmp,r_uni)))
abline(v=r_uni,lwd=2,col="green",lty=2)
legend("top",lwd=2,col="green",lty=2,legend="Uniforme",bty="n")
 
hist(dist_vmp,xlim=c(min(dist_vmp,r_agre),max(dist_vmp,r_uni)))
abline(v=r_ale,lwd=2,col="blue",lty=2)
abline(v=r_agre,lwd=2,col="red",lty=2)
abline(v=r_uni,lwd=2,col="green",lty=2)
legend("topright",lwd=2,col=c("blue","red","green"),lty=2,legend=c("Aleatório","Agregado","Uniforme"),bty="n")
 
##Observado para uma distribuição ao acaso
mean(dist_vmp)
 
##Valor Predito
E_r<-sqrt(100/10000)/2
E_r
 
#Razão observado por esperado
r_ale/E_r
r_agre/E_r
r_uni/E_r

Distância até o vizinho mais próximo.

A gente ja falou de point process aqui, agora vamos ver como é a distribuição do vizinho mais próximo, e o efeito de borda que podemos ter nele, isso é meio confuso, mas vamos la.

Suponha que nos temos o problema de desenhar linhas para unir os pares de vizinhos mais próximos em um dado conjunto de pontos (x,y) que estão mapeados em duas dimensões.

Então esse é o nosso conjunto:

01

Nos precisamos de três passos.

  • Calcular a distância para o vizinho mais próximo
  • Identificar o vizinho mais próximo de cada ponto
  • usar essas distâncias minimas para identificar todos os vizinhos mais pŕoximos

Para o primeiro passo, é só calcular a distância, usando o bom e velho triângulo retângulo.

triganguloretangulo

Então a gente implementa a fórmula d=\sqrt{(y_2-y_1)^2+(x_2-x_1)^2}

1
2
3
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}

Agora nos percorremos o pontos e calculamos um vetor de distâncias, que vamos chamar de d. Ai desse vetor a gente pega a menor distância, tirando o i porque não adianta a distância de um ponto para ele mesmo, e guardamos o índice também do ponto mais próximo.

1
2
3
4
5
6
7
for (i in 1:100) {
    for (j in 1:100){
        d[j] <- distancia(x[i],y[i],x[j],y[j])
    }
    r[i] <- min(d[-i])
    nn[i] <- which(d==r[i])
}

Agora podemos colocar no gráfico, uma linha unindo os pontos mais próximos.

02

Veja que alguns pontos tem a linha de mais próximo sobrepostos. Agora vamos identificar quem está mais próximo da borda. A borda inferior e da esquerda são zero, as bordas superiores e da direta são 1, a gente pode calcular a distância para todos os pontos, vendo o menor entre todas as bordas.

borda <- pmin(x,y,dist_borda_superior,dist_borda_direita) Veja que isso funciona, porque o valor mínimo é zero, ai temos a distancia para as bordas superiores e direta, e pegamos o valor mínimo desses conjuntos, como temos um quadrado, não tem porque pegar uma distância diagonal. Agora podemos ver quantos desses valores são menores que a distância para o vizinho mais próximo.

> sum(borda

Podemos identificar esses caras no mapa

03

Agora olha o que acontece com esses valores, quando vemos a média deles com e sem esses indivíduos da borda.

04

Comumente, uma amostra é parte de um todo, então existe mais pontos, fora dessas bordas que não vemos, se pensarmos assim, veja que alguns pontos estão longe de outros pontos meramente porque não podemos ver o ponto mais próximo dele, que está fora do nosso campo de visão, fora da borda. Então eles terem recebido uma distância maior tem do vizinho mais próximo tem haver com isso na verdade. Veja que se identificamos esses caras perto demais da borda, podemos desconfiar deles, e se retiramos eles, a gente diminui a média de distâncias do vizinho mais próximo. Mas perdemos um pouco de pontos amostrados fazendo isso.

Agora pensando numa amostragem real, esse efeito pode ser muito maior se a gente usar um retângulos como área de amostra, o famoso transecto. Mas pode ser minimizado usando um circulo, que é mais difícil de implementar no campo. Então isso é algo a se levar em conta num delineamento, quando o espaço é importante, a distribuição espacial, a localização das amostras.

Bem é isso ai, 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:
Michael J Crawley 2007 – R Book 2nd. Editora Wiley 1072pp .

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
set.seed(2)
 
##Distribuição de pontos aleatórios
x <- runif(100)
y <- runif(100)
 
##Num mapa
plot(x,y,pch=19,col="black",main="Mapa")
 
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
##Menor Distância para cada ponto
 
##Vetor para guardar a menor distância
r <- vector()
##Vetor para guardar o indice do ponto mais próximo
nn <- vector()
#Vetor auxiliar para guardar todas as distâncias
d <- vector()
 
##Para cada ponto i
for (i in 1:100) {
    ##Calcular a distancia dele até o ponto j
    for (j in 1:100){
        d[j] <- distancia(x[i],y[i],x[j],y[j])
    }
    ##Qual a menor distância? Note que tiramos o i, porque a
    ##distância de um ponto pra ele mesmo não interessa, é 0
    ##Anotamos a distância
    r[i] <- min(d[-i])
    ##E o indice do ponto correspondente
    nn[i] <- which(d==r[i])
}
 
##Desenhando uma linha até o ponto mais próximo
plot(x,y,pch=19,col="black")
for (i in 1:100){
    lines(c(x[i],x[nn[i]]),c(y[i],y[nn[i]]),col="red",lty=3,lwd=2)
}
 
##Descobrindo quem está mais próximo da borda do que do vizinho mais próximo:
 
##bordas
dist_borda_superior <- 1-y
dist_borda_direita <- 1-x
 
##Distância para a borda
borda <- pmin(x,y,dist_borda_superior,dist_borda_direita)
 
##Quantos pontos estão mais próximos da borda
##Do que do vizinho mais próximo
sum(borda<r)
 
##
plot(x,y,pch=19)
indice <- which(borda<r)
points(x[indice],y[indice],col="red",cex=1.5,lwd=2)
 
##
layout(matrix(1:4,ncol=2,nrow=2))
hist(borda,xlim=c(0,0.5),ylim=c(0,20),breaks=seq(0,0.5,length.out = 20),main="Distância para borda")
hist(borda[-indice],xlim=c(0,0.5),ylim=c(0,20),breaks=seq(0,0.5,length.out = 20),main="Sem efeito de borda")
hist(r,xlim=c(0,0.15),ylim=c(0,20),breaks=seq(0,0.15,length.out = 20),main=paste("Vizinho Mais Próximo\n Média:",round(mean(r),4)))
hist(r[-indice],xlim=c(0,0.15),ylim=c(0,20),breaks=seq(0,0.15,length.out = 20),
     main=paste("Sem efeito de borda\n Média:",round(mean(r[-indice]),4)))

Existem infinitos números primos!

Essa é a pergunta que foi respondida pelo Euclides, cerca de 350 antes de Cristo, e hoje, até hoje a gente (pelo menos eu) apanha para entender.

Vamos olhar

Teorema: Existem infinitos números primos.
Suponha que exista uma quantidade finita de números primos. Seja a lista desses números p_1,p_2,p_3,\dots,p_n a lista de todos os números primos e um m um número tal que m = p_1 p_2 p_3 \dots p_n + 1. Note que m não é divisivel por p_1 ja que dividir m por p_1 vai dar um quociente de p_2 p_3 \dots p_n e um resto de 1. Similarmente, m não é divisivel por p_2,p_3,\dots,p_n.

Agora considerando o fato que todo número maior que 1 ou é um número primo, ou pode ser escrito como um produto de números primos. É claro que m é maior que 1, já que ele é alguma coisa mais 1, então m é primo ou deve ser um produto de primos.

Suponha que m seja primo, note que m é maior que todos os números na lista p_1,p_2,p_3,\dots,p_n, então nos encontramos um número que é primo que não está na lista, contradizendo nossa suposição original de que tinha a lista de todos os números primos.

Agora suponha que por outro lado m não é primo, então ele tem que ser um produto de primos. Suponha que q é um dos primos desse produto. Então m tem que ser divisível por q. Mas a gente já viu que m não é divisível por nenhum número na lista p_1 , p_2,\dots , p_n, então denovo nos temos uma contradição com a suposição de que temos uma lista de todos os números primos.
Já que a suposição de que existem uma quantidade finita de números primos leva a uma contradição, tanto considerando m como primo ou não primo, que são as duas únicas possibilidades, então deve existir uma quantidade infinita de números primos.

Se a nossa lista de primos for 2 e 3, então 2 vezes 3 da 6, mais 1 da 7 que também é primo. Se a lista fosse 2,3 e 5, multiplicando da 30, mais 1 da 31, que também e primo, ou seja, a gente pode ir gerando primos assim. A gente não gera o próximo primo da lista, mas gera um número primo, ou seja, a sempre tem mais um número primo.

Bem é isso ai, esse negócio de prova é muito divertido, 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:
Daniel J. Velleman 2006 – How to Prove It: A Structured Approach 2ed Cambridge University Press 384pp.