Parece bem bacana, factível e útil. Acho que começar apenas com a produção das matrizes de distâncias é mais prudente (plano B). Caso consiga avançar para o modelo nulo, pense em uma função independente cuja a entrada são as matrizes geradas pela primeira função. Geralmente é o que fazemos: decompomos o problema em funções subordinadas e que depois podem ser concatenadas em uma função final. Quanto ao modelo nulo, não sei se entendi direito, mas parece que vc. quer fazer um teste de Mantel, que olha a correlação entre os valores de duas matrizes e cria um cenário nulo permutando uma delas e calculando a correlação muitas vezes (há uma solução analítica, mas gosto da permutação). Isso é legal tb, mas dá mais trabalho… posso te ajudar com isso depois que terminar a parte de gerar as matrizes!

Plano B

O plano B, mais simples, pode ser somente ler os dados, calcular as distâncias euclidianas e montar as matrizes e calcular uma correlação (Pearson ou Spearman, dependendo das características dos dados)

Cuidado para que sua função fique muito específica apenas para seus dados. Pense em como generalizar para um problema geral na área.

---

cor.esp()

Correlacao espacial entre distancias geograficas e diferencas de uma variavel de interesse

Descricao

Calcula a correlacao entre as distancias eucliadianas de um conjunto de localidades e as diferencas de uma dada variavel de interesse de localidade para localidade

Uso

cor.esp(tabela, plot.entr=TRUE, plot.hist=TRUE, metodo=“pearson”, permut=10000)

Argumentos

tabela: dataframe com 4 colunas (com titulos de colunas): coluna 1 - identificacao das localidades; colunas 2 e 3 - coordenadas geograficas x e y das localidades; coluna 4 - variavel de interesse associada as localidades (quantidade de areas nativas, riqueza de especies).

plot.entr: TRUE - retorna 2 graficos representando os dados de entrada

plot.hist: TRUE - retorna o histograma com os valores dos coeficientes de correlacao (pearson ou spearman) e o valor de p resultante das permutacoes

metodo: metodo de correlacao “pearson” ou “spearman”

permut: numero de permutacoes para o calculo do p da correlacao espacial

Detalhes

Os nomes das colunas e nomes das localidades da “tabela” serao usados para identificar titulos, eixos de graficos e dimensoes das matrizes de saida de resultados. Assim, nao e conveniente utilizar nomes grandes, para facilitar a visualizacao dos resultados.

Para a visualizacao de conjuntos de dados grandes, recomenda-se, para uma primeira visualizacao dos dados, selecionar plot.hist=FALSE e permut=100 (ou menor). Com isso a visualizacao dos dados de entrada por meio de plot.entr=TRUE sera mais rapida e apresentada maior no dispositivo grafico. Caso se deseje expandir a visualizacao dos dados de entrada, basta aumentar o tamanho do dispositivo grafico e rodar novamente a funcao.

Com plot.entr= TRUE e plot.hist=TRUE, o histograma sera apresentado juntamente com os graficos de entrada. Se quisermos visualizar somente o histograma ocupando todo o dispositivo grafico, basta selecionar plot.entr=FALSE e plot.hist=TRUE.

O histograma de saida (plot.hist) refere-se sempre ao metodo de correlacao selecionado (pearson ou spearman).

Valores

A tabela (dataframe de entrada) nao deve conter celulas sem informacao (NA) e e preciso verificar se as celulas com zeros nao representam na verdade NA. Celulas com NA ou falsos zeros podem viesar o resultado da correlacao.

O grafico Localizacao x Variavel de interesse permite visualizar a localizacao geografica das localidades (xy) e distinguir as localidades quanto aos quartis da variavel de interesse.

O grafico de correlacao espacial permite a visualizacao da correlacao entre as matrizes de distancia apresentadas no console (matriz de distancia eucliadiana entre as localidades x diferencas da variavel de interesse de localidade para localidade).

O histograma retorna as frequencias do coeficiente de correlacao (R quando Pearson e rho quando Spearman) decorrentes das permutacoes da variavel de interesse. A linha vermelha representa o valor obtido (R ou rho) e o p representa a probabilidade acumulada até o extremo da distribuicao. Quando o coeficiente (R ou rho) e negativo, o p refere-se a probabilidade acumulada entre a linha vermelha e o extremo negativo. Quando o coeficiente e positivo, p refere-se a probabilidade acumulada entre a linha vermelha e o extremo positivo. O coeficiente (R ou rho) e o p estao expressos no titulo do grafico.

Nota

Quando o metodo escolhido e “spearman”, e comum surgirem mensagens de aviso no console. As mensagens de aviso referem-se a funcao cor.test() utilizada para calcular o rho de spearman. No calculo do rho, a funcao avisa que nao e possivel retornar o p exato. Os avisos nao influenciam no calculo do rho, e o p da funcao cor.esp e obtido posteriormente por meio da permutacao dos vetores de distancia e do calculo de rho para cada uma das permutacoes.

Um numero grande de permutacoes (maior que 1000) pode fazer com que funcao demore alguns segundos para retornar os resultados. O numero padrao e de 10000 permutacoes, mas pode ser reduzido quando se deseja apenas visualizar rapidamente os graficos de entrada e matrizes de distancia.

Autor:

Alexandre Toshiro Igari

alexandre.igari@usp.br

Referencias

Eu ainda tenho que estudar um pouco mais quais o metodos que possibilitam calcular a correlacao espacial. O Ale deu a dica de estudar o metodo de Mantel, que me pareceu muito semelhante ao que eu fiz na funcao, mas preciso entender o metodo mais detalhadamente.

Veja tambem

cor.test()

Exemplos

# Criando um dataframe com dados para os exemplos:

ex.1 ← data.frame(Localidade=c(“A”, “B”, “C”, “D”, “E”), x=c(1:5), y=c(2,9,3,12,5), Riqueza=c(101,110,149,121,130))

# Funcao com matrizes de distancia (no console), todos os graficos e correlacao de pearson (no dispositivo grafico) # Para ver os graficos maiores, aumente a janela do dispositivo grafico e rode a funcao novamente:

cor.esp (ex.1, permut=200)

# So com os graficos de entrada (maiores):

cor.esp (ex.1, plot.entr=TRUE, plot.hist=FALSE, metodo=“pearson”, permut=200)

# So com o histograma de saida (maior):

cor.esp (ex.1, plot.entr=FALSE, plot.hist=TRUE, metodo=“pearson”, permut=200)

# Funcao com correlacao de spearman e 10000 permutacoes (leva alguns segundos para rodar:

cor.esp (ex.1, plot.entr=FALSE, plot.hist=TRUE, metodo=“spearman”)

---

---

#ESTRUTURA DA FUNCAO

cor.esp ← function(tabela, plot.entr=TRUE, plot.hist=TRUE, metodo=“pearson”, permut=10000)

{

nomes←colnames(tabela)

# MONTANDO AS MATRIZES DE DISTANCIA ####################################

# Matriz de distancias entre coordenadas X das localidades

dist.x ← dist(tabela[,2], method=“manhattan”,diag=TRUE, upper=FALSE)

dist.x

# Matriz de distancias entre coordenadas Y das localidades

dist.y ← dist(tabela[,3], method=“manhattan”,diag=TRUE, upper=FALSE)

dist.y

# Calculando a matriz de distancia euclidiana entre as localidades (D^2 = X^2 + Y^2)

euclid ← (dist.x^2 + dist.y^2)^0.5

dist.xy ← as.matrix(euclid)

dimnames(dist.xy) ← list(tabela[,1],tabela[,1])

dist.xy

# Calculando a matriz de distancia absoluta da variavel de interesse

dist.var ← dist(tabela[,4], method=“manhattan”,diag=TRUE, upper=FALSE)

dist.var ←as.matrix(dist.var)

dimnames(dist.var) ← list(tabela[,1],tabela[,1])

dist.var

# CALCULANDO A AUTOCORRELACAO ESPACIAL ################

# VETOR DE ORDENAMENTO DAS DISTANCIAS EUCLIDIANAS XY

# Matriz de distancias entre coordenadas X das localidades (sem zeros)

dist.x1 ← dist(tabela[,2], method=“manhattan”,diag=FALSE, upper=FALSE)

dist.x1

# # Matriz de distancias entre coordenadas Y das localidades (sem zeros)

dist.y1 ← dist(tabela[,3], method=“manhattan”,diag=FALSE, upper=FALSE)

dist.y1

# Calculando a matriz de distancia euclidiana entre as localidades (D^2 = X^2 + Y^2)

dist.xy1 ←(dist.x1^2 + dist.y1^2)^0.5

dist.xy1

# Extraindo o vetor de ordenacao das distancias euclidianas

vet.xy ← as.vector(dist.xy1)

vet.xy

length(vet.xy)

# VETOR DE ORDENAMENTO DAS DISTANCIAS COM RELACAO A VARIAVEL DE INTERESSE

# Calculando a matriz de distancia absoluta da variavel de interesse

dist.var1 ← dist(tabela[,4], method=“manhattan”,diag=FALSE, upper=FALSE)

vet.var ← as.vector (dist.var1)

vet.var

length(vet.var)

# #

if (metodo==“spearman”)

{

# CALCULO DO COEFICIENTE DE CORRELACAO DE SPEARMAN (rho)

spear←cor.test(dist.var1, dist.xy1, method=“spearman”)

spear$estimate

rho.spear←spear$estimate

rho.spear

# ALEATORIZACAO PARA O CALCULO DO p DE SPEARMAN

rho←spear$estimate

rho

for(i in 2:permut)

{

rand.var1 ← sample(dist.var1)

spear←cor.test(rand.var1, dist.xy1, method=“spearman”)

rho[i]←spear$estimate

}

if(rho.spear>=0)

{

n←rho[rho>=rho.spear]

n

length(n)

p.spear ← length(n)/permut

p.spear

}

else

{

n←rho[rho< =rho.spear]

n

length(n)

p.spear ← length(n)/permut

p.spear

}

}

# #

if (metodo==“pearson”)

{

# CALCULO DO COEFICIENTE DE CORRELACAO DE PEARSON (R)

pear←cor.test(dist.var1, dist.xy1, method=“pearson”)

pear$estimate

r.pear← pear$estimate

r.pear

# ALEATORIZACAO PARA O CALCULO DO p DE PEARSON

r←pear$estimate

r

for(i in 2:permut)

{

rand.var1 ← sample(dist.var1)

pear←cor.test(rand.var1, dist.xy1, method=“pearson”)

r[i]←pear$estimate

}

if(r.pear>=0)

{

n←r[r>=r.pear]

n

length(n)

p.pear← length(n)/permut

p.pear

}

else

{

n←r[r< =r.pear]

n

length(n)

p.pear ← length(n)/permut

p.pear

}

}

###################################################################

## CENARIO NULO

# As diferencas entre as variaveis de interesse sao aleatoriamente distribuidas pelas

# localidades

# # CENARIO ALTERNATIVO

# As diferencas entre as variaveveis de interesse estao correlacionadas com as distancias euclidianas

# entre as localidades

# # SE p<0.05, REJEITA-SE O CENARIO NULO

# Quanto maior o rho, maior a autocorrelacao espacial da variavel de

# interesse

##################################################################

#

# GRAFICOS #######################################################

#

if (plot.entr==TRUE)

{

# SEPARANDO A VARIAVEL DE INTERESSE EM QUARTIS ###########################

quartis ← quantile (tabela[,4])

quartis

str (quartis)

q25 ← quartis [2]

q50 ← quartis [3]

q75 ← quartis [4]

q25

q50

q75

pt.q1 ← tabela[tabela[,4] ⇐ q25, ]

pt.q2 ← tabela[tabela[,4]>q25 & tabela[,4]⇐q50, ]

pt.q3 ← tabela[tabela[,4]>q50 & tabela[,4]⇐q75, ]

pt.q4 ← tabela[tabela[,4]>q75, ]

pt.q1

pt.q2

pt.q3

pt.q4

#

# APRESENTANDO O GRAFICO DE COORDENADAS XY COM A VARIAVEL DE INTERESSE EM QUARTIS

x11

if (plot.hist==TRUE)

{

par(mfrow=c(2,2))

}

else

{

par(mfrow=c(1,2))

}

plot (tabela[,2], tabela[,3], bty=“l”, ylim=c(min(tabela[,3]),max(tabela[,3])+50), xlab=nomes[2], ylab=nomes[3], pch=21, main=paste(“Localizacao x ”, nomes[4]), cex.main=1 )

points (pt.q1[,2], pt.q1[,3], col=“red1”, bg=“red1”, pch=21, cex=1.2)

points (pt.q2[,2], pt.q2[,3], col=“orange1”, bg=“orange1”, pch=22)

points (pt.q3[,2], pt.q3[,3], col=“green1” , bg=“green3” , pch=23)

points (pt.q4[,2], pt.q4[,3], col=“blue1”, bg=“blue1”, pch=24)

legend (cex=0.8,(min(tabela[,2])),(max(tabela[,3])+50), bty=“n” ,title=paste(nomes[4], “ quartil”),legend=c(“0-25%”,“25-50%”, “50-75%”, “75-100%”),ncol=2, pch=c(21,22,23,24), pt.bg=c(“red1”,“orange1”,“green3”,“blue1”))

#

# APRESENTANDO GRAFICO DE DISTANCIAS EUCLIDIANAS x DIFERENCA DA VARIAVEL DE INTERESSE ENTRE AS LOCALIDADES

plot (vet.xy, vet.var, bty=“l”, xlab=“Distancia Euclidiana”, ylab=paste(“Diferenca: ”,nomes[4], names(vet.var)), main=“Correlacao Espacial”, cex.main=1)

abline(lm(vet.var~vet.xy), col=“red1”)

teste←lm(vet.var~vet.xy)

}

# APRESENTANDO O GRAFICO DO HISTOGRAMA DOS COEFICIENTES DE CORRELACAO RESULTANTES DAS PERMUTACOES

if (plot.hist==TRUE & metodo==“spearman”)

{

# Histograma Spearman

if (plot.entr!=TRUE)

{

par(mfrow=c(1,1))

}

hist (rho, ylab=“Frequencia”, xlab=“rho (Spearman)”, main=paste(“Histograma de rho\n”, “rho= ”, round(rho.spear, digits=3), “\t p= ”, round(p.spear, digits=3)), cex.main=1)

abline (v=rho.spear, col=“red1”)

}

#

if (plot.hist==TRUE & metodo==“pearson”)

{

# Histograma Pearson

if (plot.entr!=TRUE)

{

par(mfrow=c(1,1))

}

hist (r, ylab=“Frequencia”, xlab=“R (Pearson)”, main=paste(“Histograma de R\n”, “R= ”, round(r.pear, digits=3), “\t p= ”, round(p.pear, digits=3)), cex.main=1)

abline (v=r.pear, col=“red1”)

}

return(list(“Matriz de distancias euclidianas:”, dist.xy, paste(“Matriz de distancia -”, nomes[4],“:”), dist.var))

}

#

#

########################## f i m ##########################

---

---

help_funcao_cor_esp.txt cor_esp.r