O problema da percolação

Em física estatística e matemática, a teoria da percolação descreve o comportamento de clusters conectados em uma grafo aleatório. Talvez se pensarmos num exemplo mais prático, fica mais fácil entender.

Podemos dizer que a água atravessando a terra em um vaso de planta é um exemplo de percolação, a água percola através da terra. Agora lembre-se que se temos partículas muito finas e bastante aglomeradas, a água não passa, por exemplo se tivermos argila no vaso de planta, mas com a terra fofa, a água passa, afinal de conta aramos a terra não?

Mas como é a dinâmica disso? Primeiro, o vaso é um exemplo em 3 dimensões certo? Podemos, para entender melhor, representar tudo de forma mais simples, vamos pensar que a terra em um vaso 2 dimensões, como um corte no vaso.

Pense que veremos algo assim na terra:

01

Aqui, os quadrados marrons representam onde tem terra, e a parte branca é espaço vazio. Então onde está vazio, a água pode passar, pode percolar, ou seja cabe água nesses espaços brancos.

02

Cabe bastante água, temos uma terra meio fofa aqui. Primeiro vamos pensar na representação disso, veja que sabendo onde tem terra, eu consigo saber onde pode ter água, então eu preciso apenas representar onde tem terra, e nesse corte 2d, eu posso fazer isso usando uma matriz, veja que a primeira figura, pode ser representada pela seguinte matriz

> mapa [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 0 0 0 1 1 0 0 1 1 [2,] 0 1 0 0 1 1 1 0 0 0 [3,] 1 0 0 0 1 0 1 0 1 1 [4,] 0 0 0 0 1 1 1 1 0 0 [5,] 0 1 0 1 1 0 0 1 1 1 [6,] 1 0 0 1 1 1 1 1 1 1 [7,] 0 1 0 0 1 1 0 1 0 0 [8,] 0 1 0 1 1 0 0 0 0 1 [9,] 0 1 1 1 1 0 0 1 0 1 [10,] 1 0 1 1 0 1 1 1 1 0

Onde tem 1 é espaço vazio, que pode ser preenchido por água, onde tem 0 é terra.
Agora, é só pensar um pouquinho para lembrar que a água não simplesmente preenche todo o espaço. Ela desce com a gravidade, ou seja, a água é jogada na linha 1, e conforme ela tem espaço para passar, ela vai descendo para as próximas linhas.

Agora primeiro vamos tentar ver como podemos ver a percolação de forma simples, a percolação vertical.

Se a gente criar uma nova matriz de dados, do mesmo tamanho da que usamos para dizer onde tem terra (chamando ela de full)

1
full<-matrix(0,ncol=ncol(mapa),nrow=nrow(mapa))

Podemos modelar onde a água vai da seguinte forma. Primeiro a gente vai pegar todos os espaços vagos da primeira linha, copiar nessa matriz full.

1
2
3
4
##Percolação Vertical
for (j in 1:ncol(mapa)) {
    full[1,j] <- mapa[1,j]
}

O que é bem simples, depois, para todas as linhas posteriores, vamos ver se existe água acima, se tem água na linha de cima e tem espaço para ter água no espaço de baixo, preenchemos nossa matriz auxiliar com 1, e assim sabemos por onde a água vai descer.

1
2
3
4
5
for (i in 2:nrow(mapa)) {
    for (j in 1:ncol(mapa)) {
        full[i,j] <- mapa[i,j]==1 & full[i-1,j]==1
    }
}
> full [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 0 0 0 1 1 0 0 1 1 [2,] 0 0 0 0 1 1 0 0 0 0 [3,] 0 0 0 0 1 0 0 0 0 0 [4,] 0 0 0 0 1 0 0 0 0 0 [5,] 0 0 0 0 1 0 0 0 0 0 [6,] 0 0 0 0 1 0 0 0 0 0 [7,] 0 0 0 0 1 0 0 0 0 0 [8,] 0 0 0 0 1 0 0 0 0 0 [9,] 0 0 0 0 1 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 0 0

E voá-la, para saber se a água percola verticalmente, basta olhar se existe algum 1 na última linha, como não tem a água, não percola para esse mapa em questão. Claro que olhar essa matriz é meio estranho, então vamos tentar ver como aquelas primeira figuras que estávamos fazendo.

03

Veja que, se pensarmos um pouquinho, não estamos simulando direito a água, ela não vai reto, a água faz curvas, se adapta, e desse jeito ela não está indo para os lados.

Agora se pensarmos um pouquinho da para fazer isso certo? Modificando nosso primeiro código, podemos ir olhando além da linha superior, se as células adjacentes cabem água e …. Bem, já adianto que com aquele código inicial, seguindo aquele tipo de estrategia mais naive, não vamos resolver esse problema tão facilmente, porque sempre você vai ter que voltar na matriz para preencher células de antes, e vai ter que passar muitas vezes olhando as possibilidades de por onde a água pode percolar. Mas existe um jeito chamado Depth-first search. Basta pensarmos que as células da matriz onde tem água são vértices de um grafo, e estar adjacente é ter uma aresta entre esses vértices.

Ai, é so fazer uma busca profunda para cada buraco onde tem água, com a seguinte função:

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
flow<-function(open,full,i,j) {
    ##casos bases
    ##Não posso olhar fora do tamanho das linhas
    if (i<1 || i>ncol(open)){
        return(full)
    }
    ##Nem fora do tamanho das colunas
    if (j<1 || j>ncol(open)){
        return(full)
    }
    ##Se o local esta fechado, acabou a percolação
    if (open[i,j]==0){
        return(full)
    }
    ##Se eu ja mexi nessa celula, não precimo mais olhar
    if (full[i,j]==1){
        return(full)
    } ## Ja marcado como cheio
 
    ##Se eu não estou em nenhum dos casos bases, ou seja, o local não esta bloqueado
    ##e eu ainda nao observei, quer dizer que tenho que preencher com 1 aqui, percolar
    full[i,j] <- 1;
 
    ##Agora eu olho em todas as direções, se eu percolo
    ##Abaixo
    full<-flow(open, full, i+1, j);
    ##Direita
    full<-flow(open, full, i, j+1);
    ##Esquerda
    full<-flow(open, full, i, j-1);
    ##Acima
    full<-flow(open, full, i-1, j);
}

Isso é um pouco difícil de explicar sem tentar entender grafos, mas podemos olhar na figura, o que acontece.

04

Agora veja que conseguimos simular melhor o efeito da água indo para os lados, e até subindo quando é possível, estamos percolando para cima, baixo, esquerda e direita, mas veja que é simples adicionar as diagonais também, ao fazermos as chamadas recursivas. E se olharmos nossa matriz full, vemos que agora na última linha existem 1, então com esse modelo consideramos que existe uma percolação nesse mapa de terra.

Bem, como eu não expliquei a função flow direito, vamos pelo menos fazer uma animaçãozinha para ver ela funcionando que é muito legal.

animated

E é isso ai, essa é uma métrica de ecologia de paisagens, e no próximo post, vamos explorar melhor porque ela pode ser útil para descrever paisagens, 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:
Robert Sedgewick Kevin Wayne 2007 Introduction to Programming in Java: An Interdisciplinary Approach Pearson Addison Wesley 736pp

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
81
82
83
84
###
#####################################
## Problema da Percolação
#####################################
##Gerando um exemplo
set.seed(123)
p<-0.5
mapa<-matrix(sample(x=c(0,1),size=100,replace=T,prob=c(p,1-p)),ncol=10,nrow=10)
mapa
 
##Figura 1
image(t(mapa)[,ncol(mapa):1],col=c("brown","white"))
 
##Figura 2
image(t(mapa)[,ncol(mapa):1],col=c("brown","blue"))
 
 
full<-matrix(0,ncol=ncol(mapa),nrow=nrow(mapa))
 
##Percolação Vertical
for (j in 1:ncol(mapa)) {
    full[1,j] <- mapa[1,j]
}
for (i in 2:nrow(mapa)) {
    for (j in 1:ncol(mapa)) {
        full[i,j] <- mapa[i,j]==1 & full[i-1,j]==1
    }
}
 
##Figura 3
par(mfrow=c(1,2))
image(t(mapa)[,ncol(mapa):1],col=c("brown","blue"))
saida<-full
saida[!mapa==full]<-2
image(t(saida)[,ncol(mapa):1],col=c("brown","blue","white"))
 
 
### Deep Searsh
flow<-function(open,full,i,j) {
    ##casos bases
    ##Não posso olhar fora do tamanho das linhas
    if (i<1 || i>ncol(open)){
        return(full)
    }
    ##Nem fora do tamanho das colunas
    if (j<1 || j>ncol(open)){
        return(full)
    }
    ##Se o local esta fechado, acabou a percolação
    if (open[i,j]==0){
        return(full)
    }
    ##Se eu ja mexi nessa celula, não precimo mais olhar
    if (full[i,j]==1){
        return(full)
    } ## Ja marcado como cheio
 
    ##Se eu não estou em nenhum dos casos bases, ou seja, o local não esta bloqueado
    ##e eu ainda nao observei, quer dizer que tenho que preencher com 1 aqui, percolar
    full[i,j] <- 1;
 
    ##Agora eu olho em todas as direções, se eu percolo
    ##Abaixo
    full<-flow(open, full, i+1, j);
    ##Direita
    full<-flow(open, full, i, j+1);
    ##Esquerda
    full<-flow(open, full, i, j-1);
    ##Acima
    full<-flow(open, full, i-1, j);
}
 
##Usando o flow
full<-matrix(0,ncol=ncol(mapa),nrow=nrow(mapa))
for (j in 1:ncol(mapa)) {
    full<-flow(mapa, full, 1, j)
}
 
##Figura 4
par(mfrow=c(1,2))
image(t(mapa)[,ncol(mapa):1],col=c("brown","blue"))
saida<-full
saida[!mapa==full]<-2
image(t(saida)[,ncol(mapa):1],col=c("brown","blue","white"))

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)))