Mapas aleatórios e o número de componentes

Uma forma de entender o ambiente melhor, é tentar ver se apesar de construirmos mapas aleatórios, ainda sim conseguimos enxergar algum padrão nessa aleatoriedade, como o número e tamanho de manchas formados

Um jeito de produzir um mapa aleatório, é simplesmente percorrer uma matriz, e sortear preencher cada célula com 0 ou 1 com uma certa probabilidade.

Podemos fazer uma função simples de ler da seguinte forma:

1
2
3
4
5
6
7
8
9
10
11
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}

Então, pensando em construir mapas quadrados (menos número de linhas e colunas), nos geramos uma matriz de zeros de um tamanho n, então para cada coluna de cada linha (i,j) nos sorteamos um número entre 0 e 1, com distribuição uniforme, e se for menor que o nosso threshold, nos trocamos o zero por 1, senão passamos pro próximo.

Veja que podemos fazer uma imagem com essa matriz, onde quadradinhos pretos são 1 e brancos são 0.

01

Nesse caso, usamos um threshold de 0.5, mas se pegarmos um valor menor, tipo 0.2, temos menos quadradinhos pretos.

02

Agora, agora veja que isso é uma imagem dessa matriz:

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

Veja que da para ver a semelhança dos pixels pretos com os 1, em relação a localização. Até ai tudo bem. Agora o que seria uma mancha num mapa? Uma mancha seria a união de vários 1 adjacentes.

Mais ou menos assim, pegando um pixel qualquer, se houver algum pixel nas 8 posições em volta desse pixel escolhido, esses 2 pixels são u conjunto, são uma mancha.

Entao temos um pixel e 8 pixels em volta dele

1 2 3 4 x 5 6 7 8

Se em alguma posição de 1 a 8 tem outro 1, unimos os dois, como? Ue, só pensar em um pixel como um vértice e a união então é uma aresta. E assim fazemos nosso grafo.

03

Veja que nessa figura, eu dei nome para os vértices da seguinte forma:

(linha,coluna)

Ou seja, o número de vértices é igual ao número de 1 da matriz, agora, tendo o grafo é facil, basta a gente ver os componentes do grafo. Que no pacote igraph ja tem uma função implementada para fazer isso, mas poderíamos usar o algorítimo union-find como descrito no livro do Robert Sedgewick

Usando a função do igraph, teremos:

> components(grafo) $membership (8,2) (2,3) (5,3) (5,4) (6,5) (2,6) (3,7) (3,8) (5,8) (3,9) (10,9) (1,10) 1 2 3 3 3 4 4 4 5 4 6 7 $csize [1] 1 1 3 4 1 1 1 $no [1] 7

E veja que aqui conseguimos ver tanto o número de componentes formados, como o tamanho de cada componente em vértices, que é a quantidade de pixels do elemento, ou seja, o tamanho da mancha.

Certo agora estamos olhando isso para um mapinha que geramos ao acaso, e se olharmos para muitos em diferentes probabilidades de geração. Bem podemos fazer isso, com 10 repetições para cada probabilidade de geração e em média veremos isso:

04

Primeiro, podemos ver que conforme aumentamos o threshold de geração do mapa, temos menos componente, a não ser quando ele é muito baixo. Mas é so pensar, se for 0, não teremos nenhum componente, e se for 1 teremos exatamente 1 componente, cobrindo todo o mapa. Agora o meio termo é onde as coisas acontecem, rapidamente a gente começa a ter muitos componentes, até um pico la pelo 0.2 e 0.3, variável dependendo de quando fazemos a simulação, mas sempre para por aqui. Depois disso, a quantidade de componentes so diminui.

Agora, avaliando em relação ao mesmo valor de threshold, não é difícil imaginar que, quanto maior o threshold, mais arestas aparecem, mais os pixels acabam unidos, então temos manchas, em média, cada vez maiores. Mas isso demora um tempo para acontecer, tanto em probabilidades de 0.1 como 0.2, temos manchas pequenas, isso porque, a relação não é exatamente linear, so a partir de um momento, que começamos a ter um aumento de probabilidade resultando num aumento do tamanho das manchas.

Bem é isso ai, apesar dos mapas serem gerados aleatoriamente, ainda sim temos algumas certezas, o que permite alguma previsibilidade, mesmo 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 2011 Algorithms 4th ed Addison-Wesley Professional 992 pp

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#####################################
## Conectividade
#####################################
library(igraph)
 
##Função para gerar mapas quadrados aleatorios, n é o tamanho e p é a probabilidade de threshold
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}
 
##Função para transformar um mapa em um grafo, com pixels com vértices e pontos adjacentes ligados por arestas
matriz_para_grafo<-function(matriz){
    aresta<-1
    grafo<-data.frame(origem=vector(),destino=vector())
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){        
            if(matriz[i,j]!=0){            
                for(m in -1:1){
                    for(n in -1:1){
                        if(i+m>0 & i+m<=nrow(matriz) &
                           j+n>0 & j+n<=ncol(matriz)){
                            if(matriz[i+m,j+n]!=0){
                                grafo[aresta,]<-c(paste("(",i,",",j,")",sep=""),paste("(",i+m,",",j+n,")",sep=""))
                                aresta<-aresta+1
                            }
                        }
                    }
                }
            }
        }
    }
    indice<-!grafo[,1]==grafo[,2]
    indices<-which(mapa!=0,arr.ind = T)    
    grafo<-simplify(graph.data.frame(grafo[indice,],directed=F,vertices=paste("(",indices[,"row"],",",indices[,"col"],")",sep="")))  
    return(grafo)
}
 
##Função para fazer uma imagem da matriz como vemos ela na tela
plota_mapa<-function(mapa){
    image(t(mapa)[,ncol(mapa):1],col=0:1)
    }
 
 
#####################################
## Mapas aleatorios
#####################################
set.seed(123)
 
##Gera um mapa com threshold de 0.5
mapa<-gera_mapa(10,0.5)
mapa
#jpeg("01.jpg")
plota_mapa(mapa)
#dev.off()
 
##Threshold de 0.2
mapa<-gera_mapa(10,0.2)
#jpeg("02.jpg")
plota_mapa(mapa)
#dev.off()
 
##Transformando o mapa em grafo
grafo<-matriz_para_grafo(mapa)
 
#jpeg("03.jpg")
plot(grafo,vertex.size=20)
#dev.off()
 
##Usando a função do igraph
components(grafo)
 
##Realizando uma simulção com 10 repetições para probabilidades de 0.1 a 0.9
repeticoes<-10
probs<-seq(0.1,0.9,0.05)
tamanho<-15
 
##Matrizes para receber os resultados da simulação
n_compomentes<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
maior_componente<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
 
##Realizando a simulação
for(i in 1:length(probs)){
    for(j in 1:repeticoes){
        mapa<-gera_mapa(tamanho,probs[i])
        grafo<-matriz_para_grafo(mapa)
        n_compomentes[i,j]<-components(grafo)$no
        maior_componente[i,j]<-max(components(grafo)$csize)
        print(paste("Concluido i=",i," j=",j,sep=""))
    }
}
 
##Avaliando o resultado
#jpeg("04.jpg")
par(mfrow=c(2,1))
plot(probs,rowMeans(n_compomentes),type="b",xlim=c(0,1),ylab="Número médio de componentes",xlab="Threshold geração mapa",frame=F)
plot(probs,rowMeans(maior_componente),type="b",xlim=c(0,1),ylab="Tamanho do maior componente",xlab="Threshold geração mapa",frame=F)
#dev.off()

Duas espécies em metapopulações de Levins

Ja vimos metapopulações para uma espécie aqui, aqui e aqui, mas será que poderíamos adicionar mais uma espécie ao modelo de Levins?

Primeiro, vamos lembrar como é a metapopulação de levins.

\frac{dO}{dt}=cO(1-O)-eO

onde:

  • O – Ocupação
  • c – taxa de colonização
  • e – taxa de extinção

Então basicamente temos dois termos,

cN(1-N)

que faz a população crescer a uma taxa c, mas quanto maior a ocupação, mais propagulos são emitidos, mas veja que quanto maior a ocupação, menor o crescimento, que é a multiplicação por (1-N), conforme o N cresce, o valor tende a zero, porque cada vez teremos 1 menos alguma coisa, que da menos que um, estaremos multiplicando por uma fração menor.

eN

Cada mancha ocupada tem uma chance de extinção e, logo quanto mais manchas ocupadas, mais manchas veremos sendo extintas.

Agora vamos colocar outra espécie, ela vai ocupar as manchas vagas pela primeira espécie, ou seja, vagou, ela tem chance de entrar ocupar, mas a se a espécie um ocupar novamente, ela perde a competição.

Como são duas espécies, precisamos agora acompanhar a mudança na ocupação de cada uma, então precisamos de duas equações

E espécie 1 continua como a equação anterior
\frac{dO_1}{dt}=c_1O_1(1-O_1)-eO_1
mas a segunda fica
\frac{dO_2}{dt}=c_2O_2(1-O_1-O_2)-c_1*O_1*O_2-eO_2

Ou seja, a metapopulação 2 expande sua ocupação, diminuindo a expansão de acordo com o o quanto ja está ocupado pela ocupação da espécie 1 e dela mesmo, e tiramos as manchas que ela perde na competição, além do que ela perde com a extinção.

Agora, quem lembra do pacote desolve, vamos usar ele para ver como isso fica ao longo do tempo, então temos que transformar numa função para usar com a função ode.

1
2
3
4
5
6
7
8
9
10
11
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}

Agora temos que determinar por quanto tempo queremos observar a ocupação e os paramtros definir os parametros de ocupação e extinção

1
2
3
4
5
6
tempo <- seq(0, 200, by =5)
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25

E vemos o seguinte

01

A espécie um rapidamente domina o ambiente, retirando a espécie 2 da jogada.

Mesmo quando mudamos os parametros um pouco, deixando ambas com a mesma taxa de colonização, maior a taxa de extinção, lembrando que para uma espécie, quando a extinção ficava maior que a colonização, ela tendia a extinção, sendo a diferença o quão rápido isso ia acontecer.

1
2
3
4
5
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05

02

Agora se mudarmos o jogo para a espécie 1 extinguir, a espécie 2 consegue dominar, salvo que ela também não tenha a colonização menor que a extinção

1
2
3
4
5
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06

03

E existem também os casos onde podemos ter uma coexistencia

1
2
3
4
5
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05

04

Agora, as possibilidades de combinações de parâmetros são grandes, tornando mais complexo definir quando uma ou outra espécie ganha, ou quando temos a cooexistencia, mas podemos ter um feeling da coisa já. além disso, nossas equações garantem a vitoria da espécie 1, mas isso não precisa ser assim.

Bem é isso ai, lembrando um pouco de equações diferencias, que não falamos a um bom tempo aqui e vendo a parte essa expansão para a metapopulação, 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:
A. A. de Oliveira e P. I. Prado – Coexistência em Metapopulações – Roteiro no EcoVirtual
Robert May & Angela McLean 2007 – Theoretical Ecology: Principles and Applications. Oxford University 272 pp

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
##
##install.packages("deSolve")
library(deSolve)
rm(list=ls())
 
##Função
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}
 
##Tempo
tempo <- seq(0, 200, by =5)
 
##Caso 1
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
 
##jpeg("01.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 2
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("02.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 3
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("03.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 4
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05
 
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("04.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()

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

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)

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