Analise Espacial
Leonardo Castelo Branco Carvalho
23/01/2022
Partição da variância - uma espécie (adaptado de http://ecologia.ib.usp.br/)
Os três parâmetros estruturais básicos das comunidades são a riqueza, composição e abundância de espécies. Podemos tentar entender as comunidades analisando cada uma de suas partes ou podemos tentar analisar o conjunto de espécies como um todo. Sendo assim, vamos analisar, via modelagem linear, a estrutura espacial que explica a variação de uma espécie de peixe.
O conjunto de dados é composto por 3 colunas: 1) Abundância da espécie em 100 pontos de amostragem; 2) PH de cada ponto de amostragem; 3) Disponibilidade espacial de plancton para a espécie.
Carregando os dados:
dados = read.table("dados_variancia.csv", header = TRUE, dec = ",", sep = ";")
head(dados)## abundancia pH plancton
## 1 5 7.424995 0.1577539
## 2 5 7.447757 0.2513368
## 3 11 7.478305 0.3449197
## 4 10 7.510168 0.4385026
## 5 11 7.539322 0.5320855
## 6 7 7.559655 0.6256684tail(dados)## abundancia pH plancton
## 95 19 7.646624 0.5320855
## 96 20 7.694641 0.6256684
## 97 28 7.689214 0.7192513
## 98 26 7.719794 0.8128342
## 99 36 7.760034 0.9064171
## 100 21 7.800000 1.0000000A abundância de uma espécie pode variar muito, no tempo e no espaço. A pergunta fundamental aqui é o que causa essa variação. A partição de variância é um método para avaliar o quanto da variação de uma quantidade pode ser atribuída à variação de uma outra quantidade. Ou seja, a hipótese a ser testada é se a variação de distribuição da espécie em uma área varia em função da variação de um fator ambiental qualquer.
Abundância da espécie
O que causa esse padrão de distribuição?
Hipóteses:
1) Variação de PH influencia na distribuição da espécie
Neste mapa, cores claras são valores mais baixos de pH e cores escuras valores altos.
2) Variação da disponibilidade de alimento (plâncton) influencia na distribuição da espécie
Neste mapa, cores claras são regiões de baixa disponibilidade de plâncton e cores escuras alta disponibilidade.
Testando a Hipótese 1
Para avaliar a primeira hipótese vamos fazer um gráfico da relação entre a abundância da espécie e o valor do pH em cada parcela. Em seguida vamos estabelecer um modelo de regressão.
Distribuição:
plot(abundancia ~ pH, data = dados, col = "blue", pch = 15)
Regressão:
regressao.pH <- lm(abundancia ~ pH, data = dados)plot(abundancia ~ pH, data = dados, col = "blue", pch = 15)
abline(regressao.pH, col = "red", lwd = 2,)
A regressão linear nos retorna muita informação importante sobre a relação que hipotetizamos. O coeficiente de determinação \(R^2\) expressa a proporção da variação total da variável resposta (a abundância de peixes, no caso) que é “explicada” pela variação da variável explanatória pH.
Extraindo o coeficiente de determinação do modelo:
summary(regressao.pH)$r.squared## [1] 0.6857125O que mostra que 68,6% da variação da abundância da espécie é explicada por uma relação linear com o pH da água.
Testando a Hipótese 2
Para testar esta nova hipótese ajustamos a regressão da abundância em função da variável plancton
Dispersão dos dados:
plot(abundancia ~ plancton, data= dados, xlab = "Dist. espacial",
ylab= "Abundância de alimento", col = "orange", pch = 2)
Modelo de Regressão:
regressao.alimento <- lm(abundancia ~ plancton, data = dados)Gráfico de regressão:
plot(abundancia ~ plancton, data= dados, xlab = "Dist. espacial",
ylab= "Abundância de alimento", col = "orange", pch = 2, lwd = 2)
abline(regressao.alimento, lwd = 2, col = "red")
Variação explicada pelo modelo:
summary(regressao.alimento)$r.squared## [1] 0.4916512Partição da variação
Resumindo o que encontramos até aqui: uma variável ambiental (pH) explica 68,6% da variação da abundância da espécie e outra variável de estrutura espacial que explica 49,2% desta mesma variação. Como é possível que a soma dos percentuais de variação explicada passe de 100% ?
Assim, a variação explicada pelo pH contém uma parte de variação espacial na disp. de alimento. Da mesma forma, uma parte da variação explicada pela variável espacial (dist. alimento - plancton) contém efeito do pH. A essa parte da variação compartilhada chamamos variação ambiental estruturada no espaço.
Modelo Estatístico:
\(y = X + W + d\)
Medindo a relação entre as duas variáveis:
plot(pH ~ plancton, data = dados, xlab = "Dist. de alimento", col = "red", pch = 19)
Variação não explicada pelas variáveis ambientais - \(d\).
A variação não explicada por nenhuma variável é obtida a partir de uma regressão que inclui todas as variáveis (no nosso caso, pH e plancton):
regressao.geral <- lm(abundancia ~ pH + plancton, data = dados)O coeficiente de determinação corresponde à parte explicada pelas duas variáveis juntas:
var.explicada = summary(regressao.geral)$r.squared
var.explicada## [1] 0.7073219O valor do \(R^2 = 71%\). Sendo assim, a variação não explicada é obtida por:
1 - var.explicada## [1] 0.2926781Portanto, 29,3% da variação na distribuição dos peixes observada, não é explica por nenhuma das variáveis mensuradas.
Pergunta: Quais seriam os próximos passos?
Análise de autocorrelação espacial
Quando se pretende fazer algum tipo de análise com dados que são espacialmente arranjados, uma das primeiras coisas que se faz é verificar e quantificar o grau de autocorrelação espacial (ou apenas correlação espacial) da variável de interesse existente entre as unidades amostrais vizinhas.
Os correlogramas são gráficos que associam os valores dos coeficientes de autocorrelação a cada intervalo de distância nos dados e são importantes para examinar padrões de autocorrelação espacial nos dados ou nos resíduos de um modelo. Eles mostram quão correlacionados são os pares de observações espaciais quando se aumenta a distância (lag) entre eles.
O coeficiente I de Moran (Moran I) é um dos principais coeficientes de autocorrelação espacial propostos na literatura, e pode ser obtido para cada classe de distância d:
\(I(d) = [\frac{n}{W_{d}}] \frac{ \sum_{i=1}^n \sum_{j=1}^{n}w_{ij}(d)(x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^n (x_i - \bar{x})^2}\)
Em que:
\(n\) é o número de unidades; \(w_{ij}(d)\) é a matriz de conectividade da classe de distância d; \(W_{d}\) a soma de todos os \(w_{ij}(d)\), sendo o número de pares de locais por classe de distância.; \(x_i\) e \(x_j\) são os valores da variável de interessse na unidades (locais) i e j.
Portanto, o valor do índice de autocorrelação para a distância d é a média dos valores de autocorrelação espacial para cada local em toda a área de estudo. O índice varia entre -1 e 1, sendo autocorrelação positiva indicadas por valores positivos e autocorrelação negativa por valores negativos. O valor esperado para ausência de autocorrelação espacial está próximo a 0.
Os correlogramas são analisados principalmente pelos seus formatos, e estes podem dar uma ideia da estrutura espacial dos dados. De acordo com P. Legendre e Legendre (2012), alguns formatos podem ser identificados como segue:
Dessa forma, pode-se realizar a analise de autocorrelação espacial de acordo com os passos a seguir.
Carregando os dados:
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-7data(mite.xy) # coordenadas espaciais
data(mite.env) # variáveis ambientais
head(mite.xy)## x y
## 1 0.2 0.1
## 2 1.0 0.1
## 3 1.2 0.3
## 4 1.4 0.5
## 5 2.4 0.7
## 6 1.8 0.9head(mite.env)## SubsDens WatrCont Substrate Shrub Topo
## 1 39.18 350.15 Sphagn1 Few Hummock
## 2 54.99 434.81 Litter Few Hummock
## 3 46.07 371.72 Interface Few Hummock
## 4 48.19 360.50 Sphagn1 Few Hummock
## 5 23.55 204.13 Sphagn1 Few Hummock
## 6 57.32 311.55 Sphagn1 Few HummockA função correlog do pacote ncf aplica o Moran I para dados univariados, é simples e permite calcular o nível de significância dos índices por reamostragem (argumento resamp).
library(ncf)##
## Attaching package: 'ncf'## The following object is masked from 'package:vegan':
##
## mantel.correlogncf.cor <- ncf::correlog(mite.xy[,1],mite.xy[,2], mite.env$SubsDens,
increment = 0.7, resamp = 1000, quiet = T)
ncf.cor## $n
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 111 276 391 336 293 225 190 168 160 98 83 44 28 12
##
## $mean.of.class
## 1 2 3 4 5 6 7 8
## 0.4652707 1.0551781 1.7487835 2.4230810 3.1239874 3.8292778 4.5230624 5.2379728
## 9 10 11 12 13 14
## 5.9235558 6.6385055 7.3274790 7.9915320 8.6596560 9.3652742
##
## $correlation
## 1 2 3 4 5
## 0.3038757421 -0.0650338129 0.0017880755 -0.0730801720 -0.1152368673
## 6 7 8 9 10
## -0.0169742705 0.1169921938 -0.0007411494 0.0972472219 0.0617981255
## 11 12 13 14
## -0.1407385493 -0.0361292015 -0.5943652720 -0.2623354355
##
## $x.intercept
## (Intercept)
## 0.5692635
##
## $p
## [1] 0.000999001 0.200799201 0.344655345 0.122877123 0.032967033 0.490509491
## [7] 0.026973027 0.414585415 0.062937063 0.192807193 0.090909091 0.410589411
## [13] 0.003996004 0.132867133
##
## $call
## [1] "ncf::correlog(x = mite.xy[, 1], y = mite.xy[, 2], z = mite.env$SubsDens, "
## [2] " increment = 0.7, resamp = 1000, quiet = T)"
##
## attr(,"class")
## [1] "correlog"plot(ncf.cor)
abline(h = 0, lty = 2)
Correlograma de Mantel
Qual a diferença entre os índices?
Mantel é utilizado para dados multivariados!
Dessa forma pode-se calcular o correlograma de Mantel da seguinte forma:
- Carregando pacote vegan, arquivo (mite) e padronizando as variáveis:
library(vegan)data(mite)
data(mite.xy)
mite_pad <- decostand(mite, method = "max")
head(mite_pad)## Brachy PHTH HPAV RARD SSTR Protopl MEGR MPRO
## 1 0.40476190 0.625 0.13513514 0.23076923 0.3333333 0.07692308 0.2352941 1
## 2 0.04761905 0.875 0.43243243 0.00000000 1.0000000 0.00000000 0.2352941 1
## 3 0.09523810 0.375 0.02702703 0.07692308 0.3333333 0.00000000 0.1764706 0
## 4 0.54761905 0.875 0.27027027 0.15384615 0.3333333 0.00000000 0.2352941 0
## 5 0.11904762 1.000 0.35135135 0.69230769 0.0000000 1.00000000 0.0000000 0
## 6 0.45238095 0.875 0.13513514 0.69230769 0.5000000 0.15384615 0.1764706 0
## TVIE HMIN HMIN2 NPRA TVEL ONOV SUCT LCIL
## 1 0.2857143 0.02777778 0.20 0.1 0.4047619 0.05479452 0.1428571 0.069156293
## 2 0.0000000 0.00000000 0.05 0.3 0.5000000 0.36986301 0.1904762 0.190871369
## 3 0.0000000 0.00000000 0.30 0.3 0.4761905 0.23287671 0.1587302 0.123098202
## 4 0.1428571 0.05555556 0.50 0.0 0.4285714 0.64383562 0.2698413 0.149377593
## 5 0.0000000 0.08333333 0.70 0.3 0.7619048 0.58904110 0.4285714 0.006915629
## 6 0.0000000 0.55555556 0.80 0.2 0.3095238 0.52054795 0.6190476 0.004149378
## Oribatl1 Ceratoz1 PWIL Galumna1 Stgncrs2 HRUF Trhypch1 PPEL
## 1 0.17647059 0.2 0.125 1.000 0.0000000 0.0000000 0.00000000 0.0000000
## 2 0.35294118 0.0 0.125 0.375 1.0000000 0.3333333 0.03448276 0.3333333
## 3 0.17647059 0.0 0.250 0.125 0.8888889 0.0000000 0.10344828 0.0000000
## 4 0.58823529 0.2 0.000 0.125 0.2222222 0.3333333 0.06896552 0.3333333
## 5 0.05882353 0.0 0.625 0.250 0.1111111 0.0000000 0.03448276 0.0000000
## 6 0.29411765 0.0 0.125 0.125 0.8888889 0.0000000 0.13793103 0.0000000
## NCOR SLAT FSET Lepidzts Eupelops Miniglmn LRUG PLAG2 Ceratoz3
## 1 0.0000000 0.00 0.0000000 0.0000000 0 0 0 0 0
## 2 0.2857143 0.25 0.1666667 0.3333333 0 0 0 0 0
## 3 0.2857143 0.00 0.6666667 0.0000000 0 0 0 0 0
## 4 0.4285714 0.25 1.0000000 0.0000000 0 0 0 0 0
## 5 0.0000000 0.00 1.0000000 0.6666667 0 0 0 0 0
## 6 0.1428571 0.00 0.8333333 0.0000000 0 0 0 0 0
## Oppiminu Trimalc2
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0Para relacionar todas as variáveis é necessário a modelagem via modelo linear de regressão, sendo assim, deve ser realizada uma análise de Regressão entre as espécies e os dados de coordenadas dos locais de coleta:
mite.pad.lm <- lm(as.matrix(mite_pad) ~ ., data=mite.xy)Em seguida, devem ser estraídos os resíduos do modelo:
mite.pad.resid = resid(mite.pad.lm)Enfim, é calculada a matriz de distâncias:
mite.D <- dist(mite.pad.resid)De posse da matriz, é obtido o correlograma de Mantel:
mite.correlog <- vegan::mantel.correlog(mite.D, XY=mite.xy, nperm=50)Correlograma:
plot(mite.correlog)