Diagrama de Voronoi

Um diagrama de Voronoi é uma forma de partição, de partir o espaço para um conjunto de elementos. Então esse é um jeito para, por exemplo, se a gente tem um conjunto de pontos no plano, dividir esse plano deixando cada ponto dentro de um polígono adjacente aos outros pontos próximos. Ele é interessante porque tem várias aplicações, entre elas, para a epidemiologia, tendo sido usado pelo John Snow para avaliar os casos de cólera e mais coisas.

Existe muitas formas de criar o diagrama de Voronoi, mas nada melhor que começar com a forma mais naive.

Suponha que temos um conjunto de pontos assim:

1
2
x<-c(2,4,6,6)
y<-c(3,3,5,2)

Que da…

01

Agora o que a gente quer fazer, a gente quer criar um polígono que separa todos os pontos, começando com o ponto do meio, que é o segundo ponto.

Para cada ponto do conjunto, começando pelo primeiro, a gente traça uma linha entre nosso ponto que estamos trabalhando (o do meio, o segundo no vetor) e o próximo ponto.

02

Dai achamos o meio entre esses dois pontos

03

E a reta perpendicular a linha que usamos para achar o ponto, traçamos uma outra reta.

04

Agora aqui temos que tomar cuidado, porque uma linha perpendicular a uma inclinação b, numa reta a+bx, vai ser 1/b, ai aqui a inclinação da reta é zero, o que gera o problema que 1/0 não tem como resolver, então nesse caso, da para trocar o zero por um número, muito muito pequeno, tipo assim:

1
2
3
bp<- -(ifelse(is.infinite(1/b),2^64,1/b))
ap<- -mean(x[c(1,ponto)])*bp+mean(y[c(1,ponto)])
abline(ap,bp)

Certo, agora isso foi para o ponto de referência e o primeiro ponto, agora vamos repetir isso para todos os pontos

Assim calculamos a reta perpendicular05

E agora o polígono desse ponto do meio vai ser a interseção entre todas essas retas perpendiculares.

06

Agora com exceção das bordas, que vai ter um tratamento especial, basicamente é sair fazendo isso para todos os pontos. Apesar de com pouco pontos assim, ficar algo estranho.

Com muitos pontos…

07

Fica bem separadinhos. Agora é interessante observar que se a gente olhar o tamanho desses polígonos, eles vão ter propriedades tipo a que vimos aqui.

Bem é isso ai, eu comecei a olhar esse diagrama de Voronoi porque ele tem aplicações muito legais, como o algorítimo de Monmonier, mas eu nem diagrama de Voronoi sabia fazer, por enquanto pode parecer meio perdido, mas juntando com mais coisas, vai fazer mais sentido o que fizemos aqui. 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.

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
x<-c(2,4,6,6)
y<-c(3,3,5,2)
 
##01
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,10))
 
 
ponto<-2
 
##02
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,10))
lines(x[c(1,ponto)],y[c(1,ponto)],lty=3)
 
##03
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,10))
lines(x[c(1,ponto)],y[c(1,ponto)],lty=3)
points(mean(x[c(1,ponto)]),mean(y[c(1,ponto)]),pch=19,col="red")
 
##04
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,10))
lines(x[c(1,ponto)],y[c(1,ponto)],lty=3)
points(mean(x[c(1,ponto)]),mean(y[c(1,ponto)]),pch=19,col="red")
 
b<- abs(y[ponto]-y[1])/abs(x[ponto]-x[1])
a<- -b*x[ponto]+y[ponto]
abline(a,b)
 
bp<- -(ifelse(is.infinite(1/b),2^64,1/b))
ap<- -mean(x[c(1,ponto)])*bp+mean(y[c(1,ponto)])
abline(ap,bp)
 
###05
plot(x,y,pch=19,xlim=c(0,10),ylim=c(-3,10))
 
for(i in (1:length(x))[-ponto]){
 
    lines(x[c(i,ponto)],y[c(i,ponto)],lty=3)
    points(mean(x[c(i,ponto)]),mean(y[c(i,ponto)]),pch=19,col="red")
 
    b<- (y[ponto]-y[i])/(x[ponto]-x[i])
    a<- -b*x[ponto]+y[ponto]
 
    bp<- -(ifelse(is.infinite(1/b),2^64,1/b))
    ap<- -mean(x[c(i,ponto)])*bp+mean(y[c(i,ponto)])
    abline(ap,bp)
}
 
###install.packages("tripack")
library(tripack)
##06
map<-voronoi.mosaic(x,y)
plot(x,y,pch=19,col="red",xlim=c(0,10),ylim=c(-3,10))
plot.voronoi(map,pch=19,add=TRUE)
 
x<-runif(50)
y<-runif(50)
 
##07
map<-voronoi.mosaic(x,y)
plot(x,y,pch=19,col="red")
plot.voronoi(map,pch=19,add=TRUE)

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

Estatística espacial – Point process

Existem três problemas essenciais para a estatística espacial são:

1. Point Process (padrões espaciais de distribuição dos indivíduos, onde os indivíduos estão)
2. Mapas de uma variável continua de resposta (kriging)
3. Resposta espacialmente explicitas afetadas pela identidade, tamanho ou proximidade dos vizinhos.

Logo, a primeira coisa a entender é os Point process, que eu não sei dar uma tradução direito, já que nem o wikipedia tem um artigo disso em português, se tem, não esta linkado com o artigo em inglês e processos de pontos ou processos pontuais me parece que fica uma tradução estranha, não sei, mas tudo bem.

Primeiro temos que ver que estamos falando aqui de dados que consistem em de conjunto de coordenadas nos eixos x e y com alguma forma na qual foram amostradas, num plot quadrado, um transecto que é um retângulo ou talvez um círculo, se bem que eu diria que ninguém amostra círculos, mas tudo bem, a vezes as amostras são em polígonos sem uma forma definida.

A primeira pergunta mais essencial é se podemos rejeitar a hipótese de completa aleatoriedade na distribuição dos pontos, ou seja, os pontos são completamente independentes um dos outros?

02

As outras possibilidades, que são os dois opostos, que são quando temos uma distribuição regular, onde os indivíduos são mais espaçados que o acaso por algum tipo de mecanismo, como competição por exemplo que elimina os indivíduos muito próximos uns dos outros

01

E a distribuição agregada, onde os indivíduos estão amontoados juntos mais que ao acaso, como quando ocorre por um processo de reprodução com uma limitação na dispersão ou devido a heterogeneidade espacial como em locais mais propícios ao desenvolvimento de plântulas, com mais recursos ou ainda locais muito ruins para viver, sem disponibilidade de água.

03

Esses são os três padrões que vem a cabeça quando falamos de distribuição espacial. Claro que nada é tudo ou nada, existe um gradiente que vai de um padrão para o outro, e pior ainda, isso pode variar dependendo da escala que fazemos nossas observações, ou da perspectiva da espécie estudada.

Mas isso quer dizer então que uma distribuição aleatória gera dados aleatórios e uma função determinística gera os outros padrões não aleatórios?

Veja que dispersão por exemplo, é um processo que é estocástico, e ainda sim pode gerar um padrão de agregação dependente do espaço.

Podemos fazer uma simulaçãozinha aqui. Partindo do seguinte princípio. Uma semente anemocórica vai sair de um ponto, que é sua mãe, e de la vai sair com uma direção e percorrer uma distância qualquer, vamos supor que estas duas sejam variáveis aleatórias com uma distribuição normal. Basicamente o lado que a semente vai é o angulo dela com a planta mãe e a distância é uma medida de distância.

Então a gente pode fazer aqui uma função para simular isso

1
2
3
4
5
6
7
dispersa_semente <- function(r) {    
    angulo <- runif(1,0,2*pi)
    raio <- runif(1,0,r)
    x <- angulo*sin(raio)
    y <- angulo*cos(raio)
    return (c(x,y))
}

Então a gente fornece um raio máximo, a capacidade de dispersão daquela planta, sorteamos uma direção e uma distância e temos a posição que a semente vai cair. Se quisermos plotar a arvore um outro ponto que não o a coordenada (0,0), é só somar isso ao x e y.

Por exemplo, escolhemos um ponto:

1
2
centro_x <- 10
centro_y <- 10

Uma capacidade de dispersão r e o número de sementes a ser dispersadas

1
2
n <- 1000
r <- 10

E vamos colocando essas sementes no mapa.

1
2
3
4
5
6
for (i in 1:n) {
    ponto <- dispersa_semente(r)
    x <- centro_x+ponto[1]
    y <- centro_y+ponto[2]
    points(x,y,pch=19)
}

05

E olha que legal, um processo estocástico, gera uma distribuição agregada ali no centro.

Podemos plotar as distancias até o centro, como retas e podemos ver melhor o que ta acontecendo.

06

Quem dispersa mais longe, consegue ficar mais longe da semente adjacente, agora quem dispersa próximo, fica muito pertinho da outra semente, por mais que o angulo e a distância sejam variáveis aleatórias uniformes.

Se a gente fizer essa simulação, e salvar os valores de angulo e raio usando para gerar os dados, e depois as posições x y e distancias da planta mãe.

07

Veja que a distancia da planta mãe ainda é uma variável aleatória uniforme, mas tanto nos eixo x quanto no y, temos uma acumulação no centro, perto da planta mãe. Ou seja um processo estocástico, aleatória que gera um padrão de distribuição não aleatório.

Acho que isso é meio confuso, o fato de um evento totalmente aleatório gerar um padrão de distribuição, mas a vida sempre é mais complexa que esperamos. E isso serve de exemplo também para pensarmos que não é porque o padrão de distribuição não é aleatório, que o processo que gera ele não possa ser aleatório.

Bem é isso ai, os posts andam meio lento, mas espero votar a estudar, ops, postar com mais frequência. 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:

M. J. Crawley 2012 – R book 1076pp Wiley

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
#####################################
## Tipos de Point Process
#####################################
 
x_coords<-c(0.3,0.6,0.8,0.2,0.4,0.6,0.8,0.3,0.6,0.8,0.2,0.4,0.6,0.8)
y_coords<-c(0.3,0.3,0.3,0.5,0.5,0.5,0.5,0.7,0.7,0.7,0.9,0.9,0.9,0.9)
jpeg("01.jpg")
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Espaçado",bty="n")
dev.off()
 
x_coords<-c(0.2,0.9,0.92,0.93,0.78,0.19,0.21,0.5,0.55,0.61,0.94,0.96,0.37,0.41,0.44)
y_coords<-c(0.1,0.2,0.22,0.23,0.3,0.36,0.37,0.4,0.39,0.39,0.43,0.44,0.56,0.53,0.56)*1.8
jpeg("02.jpg")
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Aleatorio",bty="n")
dev.off()
 
x_coords<-c(0.2,0.22,0.23,0.22,0.2,0.56,0.57,0.87,0.89,0.89,0.89,0.9,0.92,0.925,0.926)
y_coords<-c(0.2,0.22,0.22,0.225,0.28,0.7,0.71,0.8,0.81,0.82,0.83,0.83,0.83,0.82,0.81)
jpeg("03.jpg")
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Agregado",bty="n")
dev.off()
 
jpeg("04.jpg")
par(mfrow=c(3,1))
x_coords<-c(0.3,0.6,0.8,0.2,0.4,0.6,0.8,0.3,0.6,0.8,0.2,0.4,0.6,0.8)
y_coords<-c(0.3,0.3,0.3,0.5,0.5,0.5,0.5,0.7,0.7,0.7,0.9,0.9,0.9,0.9)
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Espaçado",bty="n")
x_coords<-c(0.2,0.9,0.92,0.93,0.78,0.19,0.21,0.5,0.55,0.61,0.94,0.96,0.37,0.41,0.44)
y_coords<-c(0.1,0.2,0.22,0.23,0.3,0.36,0.37,0.4,0.39,0.39,0.43,0.44,0.56,0.53,0.56)*1.8
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Aleatorio",bty="n")
x_coords<-c(0.2,0.22,0.23,0.22,0.2,0.56,0.57,0.87,0.89,0.89,0.89,0.9,0.92,0.925,0.926)
y_coords<-c(0.2,0.22,0.22,0.225,0.28,0.7,0.71,0.8,0.81,0.82,0.83,0.83,0.83,0.82,0.81)
plot(x_coords,y_coords,xlim=c(0,1),ylim=c(0,1),pch=19,main="Agregado",bty="n")
dev.off()
 
#####################################
## Pontos dentro de um circulo
######################################
 
dispersa_semente <- function(r) {    
    angulo <- runif(1,0,2*pi)
    raio <- runif(1,0,r)
    x <- angulo*sin(raio)
    y <- angulo*cos(raio)
    return (c(x,y))
}
 
##Pontos circulares
centro_x <- 10
centro_y <- 10
 
jpeg("05.jpg")
plot(centro_x,centro_y,ylab="",xlab="",ylim=c(0,2*centro_y),xlim=c(0,2*centro_x),type="n")
 
n <- 1000
r <- 10
 
for (i in 1:n) {
    ponto <- dispersa_semente(r)
    x <- centro_x+ponto[1]
    y <- centro_y+ponto[2]
    points(x,y,pch=19)
}
dev.off()
 
##retas
centro_x <- 10
centro_y <- 10
 
jpeg("06.jpg")
plot(centro_x,centro_y,ylab="",xlab="",ylim=c(0,2*centro_y),xlim=c(0,2*centro_x),type="n")
 
 
n <- 1000
r <- 10
 
x<-c()
y<-c()
a<-c()
raio<-c()
distancias<-vector()
 
for (i in 1:n) {    
    a[i] <- runif(1,0,2*pi)    
    raio[i] <- runif(1,0,r)    
    xi <- a[i]*sin(raio[i])
    yi <- a[i]*cos(raio[i])
    ponto <-c(xi,yi)
 
    x[i] <- centro_x+ponto[1]
    y[i] <- centro_y+ponto[2]
    lines(c(centro_x,x[i]),c(centro_y,y[i]))
    distancias[i]<-dist(data.frame(origem=c(centro_x,x[i]),destino=c(centro_y, y[i])))
}
dev.off()
 
jpeg("07.jpg")
layout(matrix(c(1,1,1,2,2,2,3,3,4,4,5,5),ncol=6,nrow=2,byrow=T))
hist(a,main="Raios")
hist(raio,main="Angulos")
hist(distancias,main="Distancias")
hist(x,main="Coordenadas em x")
hist(y,main="Coordenadas em y")
dev.off()

Estatística Espacial e os dados de John Snow

Sempre em biologia é recorrente algumas preocupações quanto as análises de dados. Uma delas é a correlação filogenética que ja vimos aqui.
Outra preocupação é dependência espacial dos dados. Eu acho isso bem complicado de lidar, mas ainda bem que muitas pessoas já tiveram essa preocupação e inventaram soluções bem inteligentes para lidar com a dependência espacial dos dados, que aliás não deve ser visto como um problema, mas uma rica fonte de informação.

Posso estar falando bobeira, mas me parece que quem começou com isso foi um cara chamado John Snow, que não é o Jon Snow do guerra dos tronos, mas outro Snow, o epidemiologista.

Recentemente eu vi os dados dele disponíveis em um post de um blog ai resolvi baixar para dar uma olhada, é so baixar la também.

Os dados são bem legais, as coordenadas são quadras de um bairro, o bairro de Soho em Londres. Para cada quadra, ele registrou o número de casos de ocorrência de cólera porque estava rolando uma grande epidemia na época e ele queria entender melhor como tudo estava acontecendo para poder tomar medidas eficientes de controle contra a cólera, tendo em mente que na época, as pessoas não sabiam como se contraia a cólera.

Ai qual a relação disso com estatística espacial? Bem vamos começar dando uma olhada nos dados do senhor Snow.

01

Aqui cada ponto preto representa um local onde teve casos de cólera, e o tamanho do ponto preto representa o número de casos nesse local. Os pontos azuis são onde haviam poços comunitários de aguá. Naquela época não tinha aguá encanada tratada, dai as pessoas tinha que se deslocar das suas casas, até esses poços para pegar aguá.

Olha so parece que existe uma maior ocorrência de cólera no meio ali perto de um poço. Hoje isso faz sentido fácil ja que a gente sabe como se pega cólera, mas lembre-se que na época não se sabia.

Mas continuando, uma coisa que a gente poderia fazer é ver a densidade de pontos de acordo com os locais do gráfico, isso é como fazer um histograma, so que em duas dimensões. Então organizamos os dados para proceder com essa idéia e fazemos mais um gráfico, e agora vamos dar nomezinhos ao poços.

02

Nessa figura, nos representamos a densidade com esses contornos. Eu usei um gradiente do cinza para o vermelho para representar a intensidade de casos, quanto mais vermelho, significa que o número de casos é maior, quanto mais cinza significa que o número de casos é menor. Isso ajuda a evidencia o padrão que o primeiro gráfico estava sugerindo, que existe uma incidência maior de casos numa parte do bairro, e que coincidentemente está em volta de um poço especifico, que chamamos aqui de poço 1, como apresentado nessa mesma figura.
Olha ae que legal, então tem uma região que tem mais casos, e ela é em volta de um poço.

Agora vamos fazer o seguinte, vamos pegar o ponto de cada caso de cólera, e medir a distancia dele (distancia euclidiana, uma linha reta) até os poços e veremos se da para ver alguma coisa.
Bem uma coisa que não ficou representado na figura acima, é que existem 8 poços na região, mas como alguns estão meio longe dos casos de cólera, ele nem ficaram na área do gráfico, mas temos as coordenadas deles.

Então calculamos a distância de cada poço até cada caso de cólera e vamos fazer um histograma com esses dados para cada poço e ver se existe alguma coisa acontecendo.

03

Opa olha ai que figura interessante. A primeira coisa a se notar é que o poço 1 tem a distribuição de distancias desviada para o 0.
A menor distância possível é zero certo correto? E vamos pensar, os casos estão em média com valores de distancias baixo desse poço, compare com os outros poços, parece que eles são bem mais distantes dos casos. Isso quer dizer que os doentes estão em média pertinho do poço 1. O que acontece é que o poço um tinha a água contaminada, e todo mundo ia buscar água no poço mais próximo, para andar menos uai, então quanto mais próximo, até uma certa distancia, as pessoas so pegavam água contaminada no poço 1 e ficavam doente, conforme mais distante, o cara poderia ir no poço 1 ou em outro poço com distancia igual, logo menos gente ficava doente, agora quem tava muito longe do poço 1 ia pegar água em outro poço mais próximo, logo se o outro poço não tinha água contaminada eles não ficavam doente. Isso gera esse padrão espacial onde a gente ve uma distribuição maior de casos numa área e menor em outra. E gera um correlação dos dados com alguma coisa no espaço, que é o poço contaminado.

Mas olhando somente os dados de caso de cólera, a gente ve que os dados tem uma auto-correlação. Mas que diabos é uma auto-correlação?

É quando a existe uma correlação dos dados com eles mesmo. Como uma população de passarinho crescendo, o tamanho da população na próxima medida sempre depende da medida anterior, so lembrar aqui e aqui. Nesse caso, qualquer parzinho de pontos perto do poço vão ter um valor alto de casos de cólera, qualquer parzinho de pontos longe do poço vão ter um número baixo casos de cólera. Ou seja vão ter valores próximos. Agora se a gente pegar um ponto perto do poço e um ponto longe do poço, os valores de casos de cólera deles vão ser muito diferentes uns dos outros. Isso é um tipo de dependência espacial, quanto mais próximo os pontos, mais iguais eles são, quanto mais longe um ponto do outro, mais diferente eles são, e isso é perto e longe espacialmente, de distancia mesmo que a gente tem que caminhar.

Uma maneira de representar isso graficamente é com um variograma. Mas o que é um variograma? É como é essa relação acima que eu falei.
A gente pode usar o pacote geoR para fazer isso. Aliás um dos autores desse pacote é um Brasileiro, o Paulo Justiniano, que vira e mexe esta ajudando as pessoas na lista de R Brasileira R-Br.

04

Olha ai nosso variograma, nele a gente ve, quanto mais próximos os pontos, maior a similaridade quanto ao casos de cóleras deles e quanto mais longe, menor a similaridade. Até uma hora que não tem mais relação, lembrando que essa é a relação que a gente ve aqui, ela não precisa ser somente desse jeito, mas essa figura chamada de variograma, e me perdoem se estiver usando o nome incorreto, ajuda a avaliar se existe correlações espaciais nos dados e como elas são.

Nesse caso eu fiquei pensando por um tempo se estava interpretando corretamente esse gráfico, ai resolvi fazer um teste, uma pequena simulação num padrão agregado de pontos num local e distribuídos aleatoriamente no espaço e ver como saiam as figuras.

05

A primeira coluna são dados com um padrão espacial agregado num local, a segunda coluna são dados com uma distribuição aleatória, eu gerei eles assim. Dai eu fiz o gráfico de contorno para representa a densidade deles no espaço e por ultimo o variograma para ver como ficava. Com o padrão espacial agregado ficou exatamente como os com o dado do Snow, uma dependência onde quanto mais próximo, maior a auto-correlação espacial, agora nos dados sem padrão nenhum, como deveria ser, não vemos autocorrelação nenhuma, uma linha reta.

Eu não entendo quase nada do tópico, mas eu estou tentando ler o livro do Peter Diggle e do Paulo Justiniano Ribeiro, o Model-Based Geostatistics, que é bem legal. Uma excelente leitura, que ajuda a entender como detectar padrões espaciais e levar isso em conta em modelos gerais lineares e cia. Eles são autores também do pacote geoR, que garante que tudo que a gente aprender no livro a gente vai conseguir usar, implementando com ajuda do R.

Mas eu acho legal pelo menos começar a pensar nisso, em padrões espaciais, ter uma idéia do que se trata e saber para onde correr caso surja a necessidade.

E para finalizar segue o mapinha original do John Snow, olha que legal ele fazia uns pequenos histogramaszinhos por quadra, para detectar o padrão espacial dos dados, e gente reclama que tudo é difícil hoje, imagina com lápis e papel, ai sim devia ser tenso conseguir fazer essas coisas 🙂

Mapa original do John Snow

library(MASS)
library(maptools)
library(geoR)
 
#Abrindos os dados do John Snow
pumps<-readShapePoints("Pumps.shp")
cholera<-readShapePoints("Cholera_Deaths.shp")
 
#Olhando os dados em uma figura
plot(cholera,pch=19,cex=log(cholera@data[,2]+1))
points(pumps,pch=19,col="blue")
axis(1)
axis(2)
legend("topleft",pch=19,col=c("black","blue"),
       legend=c("Casos de Colera","Poços"),bty="n")
 
#Repetindo os pontos, pelo número de casos
lat<-vector()
for(i in 1:length(cholera@data[,2])) {
  lat<-c(lat,rep(cholera@coords[i,1],cholera@data[i,2]))
}
long<-vector()
for(i in 1:length(cholera@data[,2])) {
  long<-c(long,rep(cholera@coords[i,2],cholera@data[i,2]))
}
 
#Calculando densidade de pontos
contorno<-kde2d(lat,long)
 
#Observando a densidade de casos na area
cores<-colorRampPalette(colors=c("gray","red"))
contour(contorno,col=cores(20),lwd=2,drawlabels =F,nlevels = 20,frame=F)
points(pumps,pch=19,col="blue",cex=2)
legend("topleft",pch=19,col=c("blue"),legend=c("Poços"),bty="n")
text(x=pumps@coords[1:5,1],y=pumps@coords[1:5,2]+16,1:5)
 
 
#Lista de distancias entre Poços
distancias<-list()
for(i in 1:length(pumps)) {
  distancias[[i]]<-vector()
}
#Lista de distancias entre locais de ocorrência de colera
for(j in 1:length(pumps)){
  for(i in 1:length(cholera@data[,2])) {
    distancias[[j]][i]<-dist(rbind(pumps@coords[j,1:2],cholera@coords[i,1:2]))
  }
}
 
#Lista de distancias entre poços e ocorrencias de colera
distancias2<-list()
for(i in 1:length(pumps)) {
  distancias2[[i]]<-vector()
}
for(j in 1:length(pumps)) {
  for(i in 1:length(cholera@data[,2])) {
    distancias2[[j]]<-c(distancias2[[j]],rep(distancias[[j]][i],cholera@data[i,2]))
  }
}
 
#Histogramas das distancias dos poços para ocorrências de colera
layout(matrix(1:8,byrow=T,ncol=2))
for(i in 1:8) {
  hist(distancias2[[i]],breaks=seq(0,700,by=25),ylim=c(0,100),
       main=paste("Poço ",i),xlab="Distancia do Poço",ylab="Número de casos")
}
 
 
#Variograma
breaks<-seq(0,500, by=20)
variograma<-variog(coords=cholera@coords,data=cholera@data[,2],breaks=breaks)
var.summary <- cbind(c(1:(length(breaks)-1)), variograma$v, variograma$n)
colnames(var.summary) <- c("lag", "semi-variance", "# de pares")
#Verificando número de casos por distancia
var.summary
 
#Plot do Variograma
plot(variograma,type = "b",frame=F,xlab="Distância")
 
 
#simulando alguns dados para entender melhor o variograma
set.seed(51)
#dados com dependencia espacial
xpois<-rpois(200,10)
ypois<-rpois(200,10)
tabela<-table(xpois,ypois)
dados<-data.frame(expand.grid(rownames(tabela),colnames(tabela)),
                  c(tabela))
dados<-dados[dados[,3]>0,]
dados[,1]<-as.numeric(as.character(dados[,1]))
dados[,2]<-as.numeric(as.character(dados[,2]))
 
#dados sem dependencia espacial
xr<-runif(200,min(xpois),max(xpois))
yr<-runif(200,min(ypois),max(ypois))
tabelar<-table(xr,yr)
dadosr<-data.frame(expand.grid(rownames(tabelar),colnames(tabelar)),
                   c(tabelar))
dadosr<-dadosr[dadosr[,3]>0,]
dadosr[,1]<-as.numeric(as.character(dadosr[,1]))
dadosr[,2]<-as.numeric(as.character(dadosr[,2]))
 
#Grafico com simulações
layout(matrix(1:6,ncol=2,nrow=3,byrow=F))
#Dependencia espacial
plot(x=dados[,1],y=dados[,2],cex=log(dados[,3]+1),pch=19,
     main="Com Padrão espacial",xlab="",ylab="",frame=F)
contour(kde2d(xpois,ypois),col=cores(20),lwd=2,drawlabels =F,nlevels=20,frame=F)
plot(variog(coords=dados[,1:2],data=dados[,3]),frame=F,
     ylim=c(-1,3),xlim=c(0,20))
 
#Sem dependencia espacial
plot(x=dadosr[,1],y=dadosr[,2],cex=log(dadosr[,3]+1),pch=19,
     main="Sem padrão espacial",xlab="",ylab="",frame=F)
contour(kde2d(xr,yr),col=cores(20),lwd=2,drawlabels =F,nlevels=20,frame=F)
plot(variog(coords=dadosr[,1:2],data=dadosr[,3]),frame=F,
     ylim=c(-1,3),xlim=c(0,20))