Variabilidade espacial do pH do solo da região de Unaı́-MG
Ana Julia Righetto 1
Luiz Ricardo Nakamura 1
Diogo Eberhardt 2
Paulo Justiniano Ribeiro Junior 3
1
Introdução
Segundo Druck et al. (2004), compreender a distribuição espacial de dados provenientes de
fenômenos ocorridos no espaço constitui um grande desafio para esclarecer questões centrais
em diversas áreas do conhecimento, como saúde, ambiente, geologia, agronomia, entre outras.
Tais estudos tornam-se cada vez mais comuns devido a disponibilidade de sistemas relacionados
à informação espacial, como o Sistema de Informação Geográfica (SIG) e o Sistema de Posicionamento Global (GPS, do inglês “Global Positioning System”), bem como os avanços nos
procedimentos computacionais que permitem e facilitam a análise ou representação do espaço
e dos fenômenos que nele ocorrem. Assim, é crescente o interesse em analisar e ajustar modelos para dados georrefereciados, tornando a estatı́stica espacial uma área do conhecimento em
constante ascenção.
A estatı́stica espacial tem como intuito mensurar propriedades e relacionamentos, levando
em consideração a localização espacial do fenômeno estudado de forma explı́cita. Esse ramo
da estatı́stica é dividido em três áreas de estudo: geoestatı́stica, dados de área e processos pontuais. Como o que distingue cada uma dessas categorias é o tipo de dado aleatório utilizado, é
normal que em cada uma destas três áreas, existam diferentes métodos estatı́sticos para analisar
e descrever um determinado conjunto de dados.
Neste trabalho, buscou-se avaliar a variabilidade espacial do pH da água, em um região
localizada em Unaı́, noroeste de Minas Gerais, Brasil. Este tipo de variável (contı́nua), do
ponto de vista da estatı́stica espacial, requer a utilização da geoestatı́stica, que é um ramo da
geografia matemática e da estatı́stica, que une o conceito de variáveis aleatórias com o conceito
de variáveis regionalizadas, gerando um novo conceito de funções aleatórias. Neste ramo, as
técnicas que mais se destacam são krigagem e a simulação estocástica e, com o auxı́lio destas
e outras técnicas, é possı́vel calcular um valor de uma dada região em estudo, para cada centro
da célula de uma malha tridimensional, valor este condicionado aos dados amostrados e a uma
função de correlação espacial entre estes dados. Resumidamente, Matheron (1963) define a
1 ESALQ/USP.
e-mail: [email protected]
2 DEP
3 UFPR
1
geoestatı́stica como a aplicação do formalismo de funções aleatórias para o reconhecimento e
estimação de fenômenos naturais.
Além da aplicação da geoestatı́stica em sua abordagem frequentista, neste trabalho foi utilizada a abordagem bayesiana que permite a incorporação dos conhecimentos a priori de um
pesquisador da área em estudo. Após os ajustes frequentista e bayesiano realizou-se uma
comparação entre os dois.
2
2.1
Material e métodos
Área de estudo
O estudo foi conduzido no Cerrado (savanas neotropicais brasileiros), como parte da investigação sobre o impacto das práticas de agricultura de conservação da fertilidade do solo e produtividade das culturas. A área em estudo está localizada em Unaı́, noroeste de Minas Gerais,
Brasil (Figura 1). O experimento de campo foi localizado na Faculdade Juvêncio Ferreira Martins Agrı́cola de Unaı́ (16o 32’26”S, 46o 50’44”W e altitude de 600 m). Essa região é caracterizada pelo clima subúmido tipicamente tropical do Cerrado brasileiro. Segundo a classificação
de Köppen, o clima é tropical úmido e seco ”Aw”(ou clima de savana).
A área em estudo é um quadrado de 80 m × 80 m. As amostras de solo foram recolhidas,
em 2010, por meio de uma grade de 5 × 5 m, a uma profundidade de 0-20 cm, resultando num
total de 240 amostras. Das variáveis coletadas, neste trabalho, analisou-se o pH da água.
Figura 1: Região de estudo
2
2.2
Abordagem frequentista
Quando a variável de interesse segue distribuição gaussiana, torna-se possı́vel a construção
de um variograma a partir do cálculo da semivariância γ(h) dada por:
γ(h) =
N
1
[y(xi + h) − y(xi )]2 .
∑
2N(h) i=1
(1)
Para analisar como os valores dos dados se comportam conforme várias distâncias estipuladas e assim, analisar o grau de dependência espacial da variável e definir os parâmetros
necessários para estimar caracterı́sticas em locais não amostrados, utiliza-se um gráfico de γ(h)
por h denominado de semivariograma, principal ferramenta para o diagnóstico da existência de
correlação entre os pontos em estudo.
O semivariograma analisa o grau de dependência espacial entre amostras dentro de um
campo experimental, além de definir parâmetros necessários para a estimativa de valores para
locais não amostrados, através da técnica de krigagem (SALVIANO, 1996). Ainda, segundo Ribeiro Jr. (1995), os semivariogramas são preferidos para caracterizar a estrutura de continuidade
espacial da caracterı́stica avaliada, pois exigem hipóteses de estacionaridade menos restritivas
(hipótese intrı́nseca).
O ajuste de um modelo matemático aos dados no gráfico de γ(h) por h, isto é, à uma função,
é um dos possı́veis métodos para estimar os parâmetros do semivariograma que consiste em
construir classes de pares de pontos equidistantes e realizar uma espécie de estudo da evolução
da correlação entre os grupos conforme se caminha na região em estudo, ou seja, uma função
matemática dada pela equação 1 que depende de h, a distância entre pontos de um mesmo
grupo, que fornece o efeito de sinal de σ2 , o qual mede o quanto o padrão espacial está nı́tido
- quanto maior o sinal maior a confiança na existência de dependência espacial -, além do
efeito do ruı́do τ2 , também chamado efeito pepita, que define a variabilidade não explicada.
Para afirmar a dependência espacial dos dados com mais segurança, constrói-se um envelope
de variogramas que funciona como uma espécie de intervalo de credibilidade para a aparente
dependência espacial do variograma empı́rico obtido.
Segundo Vieira (2000), em um semivariograma ajustado, o valor da semivariância aumenta
na medida em que aumenta-se a distância de separação entre os pontos, até que seja atingido
um patamar no qual ele se estabiliza. Este patamar, atingido quando a variância dos dados
fica constante em relação a distância entre os pontos amostrados, permite a determinação da
distância limite entre dependência e independência entre as amostras.
Foram ajustados dois modelos, um considerando como função de correlação a função exponencial e o outro considerando a função Matérn, que são representados respectivamente por (2)
e (3):
3
−h
γ(h) = exp
φ
γ(h) =
k h
h
Kk
(k−1)
φ
φ
2
Γ(k)
1
(2)
(3)
no qual Γ(.) é a função gama, Kk é a função Bessel de ordem k, ||h|| é a distância euclidiana
entre duas localizações e φ é o parâmetro de dependência espacial.
Outro critério utilizado para afirmar essa dependência espacial é denominado Cálculo da
Razão de Dependência Espacial (RD). O valor do RD é dado em percentagem segundo uma
função que depende do efeito pepita do modelo (τ2 ) e do sinal de σ2 (CAMBARDELLA et al.,
1994). A equação para o cálculo dessa razão é dada por:
RD =
τ2
×100%.
τ2 + σ2
(4)
A interpretação dos valores obtidos na equação (4) é realizada da seguinte maneira:
• RD ≤ 25%, indica dependência espacial forte;
• 25% < RD < 75%, a dependência é moderada;
• RD ≥ 75%, a dependência é fraca.
Os modelos utilizados foram ajustados por meio do método da máxima verossimilhança e
como critério para a escolha do melhor ajustado, foram utilizados os valores do RD, o critério de
Akaike (AIC) e o critério de Schwarz (BIC), dados pelas equações (4), (5) e (6) respectivamente.
AIC = −2ln(L) + 2p
(5)
BIC = −2ln(L) + ln(np)
(6)
Nas equações (5) e (6), L é a função de verossimilhança e n é o número de parâmetros
ajustados. Para ambos os critérios dessas equações, o modelo melhor ajustado é aquele com o
menor valor obtido de AIC e BIC.
Assim, analisando a existência da dependência espacial, pode-se estimar pontos desconhecidos na área em estudo, pois estes serão semelhantes aos seus vizinhos, aos pontos amostrados.
Esta estimação é realizada por meio da krigagem, que estima valores com variância mı́nima.
Segundo Landim (2006), krigagem é um processo de estimativa de valores de variáveis distribuı́das no espaço e/ou no tempo, a partir de valores adjacentes enquanto considerados como
interdependentes pelo semivariograma. O uso desta técnica permite a construção de uma média
4
ponderada que atribui pesos aos pontos vizinhos que são conhecidos, sendo que, usualmente,
para vizinhos mais próximos, o peso atribuido é maior. Esta interpolação de valores é dada
conforme a seguinte equação:
N
Z(x0 ) = ∑ λi Z(xi )
(7)
i=1
As construções de mapas com os valores obtidos por meio da krigagem são importantes
para a verificação e interpretação da variabilidade espacial dos dados em estudo.
2.3
Abordagem bayesiana
Em geoestatı́stica, os parâmetros dos modelos são estimados e utilizados na obtenção das
predições espaciais das variáveis em estudo. Desta forma, pode-se dizer que a incerteza sobre estes parâmetros pode ser incorporada na incerteza associada às predições realizadas. No
enfoque bayesiano a incerteza sobre os parâmetros é representada por distribuições de probabilidades. Ainda, as associações espaciais são melhores capturadas quando se faz uso de modelos
hierárquicos que constróem dependência em diferentes estágios, o que segue o paradigma bayesiano de inferência estatı́stica.
Quando se utiliza a metodologia da geoestatı́stica clássica na estimação de parâmetros,
por meio do método da máxima verossimilhança, utiliza-se apenas as informações provenientes dos dados, do modelo adotado e das teorias assintóticas. Entretanto, quando se utiliza os
métodos bayesianos, além de todos esses conhecimentos, também se utiliza informações a priori sobre os parâmetros, isto é, pode existir um conhecimento inicial que é expresso por meio
das distribuições a priori e, desta forma, na prática, não existe diferença entre parâmetros e
observações, sendo que ambos são tratados como quantidades aleatórias (RIGHETTO, 2013).
Neste trabalho, foram realizadas 10000 iterações e utilizou-se prioris não informativas para os
parâmetros, em que assumiu-se a uniforme para o parâmetro τ2 , recı́proca para os parâmetros
σ2 e φ e a normal para o parâmetro β0 .
Através do Teorema de Bayes obtém-se a distribuição a posteriori dos parâmetros, dada por:
π(θ|y) = Z
L(y|θ)π(θ)
(8)
L(y|θ)π(θ)dθ
Θ
em que θ é o parâmetro, ou vetor de parâmetros, de interesse; o denominador da função é
constante sob a distribuição de probabilidade de θ|y; L(y|θ) é a função de verossimilhança de
y|θ e π(θ) é a distribuição a priori de θ, que representa o conhecimento prévio da distribuição
de probabilidade dos parâmetros.
Segundo Kolman (2004), retirando o termo constante de (8) tem-se:
π(θ|y) ∝ L(y|θ)π(θ)
5
(9)
sendo que π(θ|y) dado em (9) não é uma distribuição de probabilidades, porém ωπ(θ|y) é,
se houver uma escolha adequada para a constante ω, que é obtida por métodos analı́ticos ou
numéricos.
Como já mencionado, após observar Y = y, existe o interesse em realizar previsões para
uma quantidade Yo não observada que, relacionada a θ é dado por:
π(yo |y) =
Z
Z
π(yo |θ, y)π(θ|y)dθ,
π(yo , θ|y)dθ =
Θ
Θ
sendo que na maioria das vezes, essa integral não possui solução analı́tica e faz-se necessário
métodos de aproximação para resolvê-la.
Depois de definir a distribuição a posteriori dos parâmetros é possı́vel resumir informações
sobre os mesmos por meio de técnicas inferenciais, utilizando, usualmente, probabilidades e
esperanças. Nesta etapa, muitas vezes é necessário utilizar métodos computacionais intensivos
como o Monte Carlo via Cadeias de Markov (MCMC). Na geoestatı́stica, o amostrador de Gibbs
e o algoritmo Metropolis-Hastings são os métodos mais utilizados.
Neste trabalho, utilizou-se modelos gaussianos geoestatı́sticos possibilitando a obtenção de
intervalos de credibilidade para os parâmetros e, com isso, foi possı́vel realizar comparações
com os resultados obtidos por meio da metodologia frequentista.
3
Resultados e discussões
Todas as análises foram realizadas no software R, utilizando o pacote geoR (RIBEIRO
JÚNIOR e DIGLLE, 2001).
Primeiramente, foi realizado o teste de Shapiro-Wilk na variável pH da água e constatou-se
normalidade ao nı́vel de 5% de significância. Na Tabela 1 são apresentadas estatı́sticas descritivas desta variável.
Tabela 1: Estatı́stica descritiva da variável pH da água
Mı́nimo 1o Quartil Mediana Média 3o Quartil Máximo
4,150
4,898
5,070
5,059
5,222
6,030
O gráfico superior a esquerda, apresentado na Figura 2, categoriza os dados nos quartis
amostrais das observações, em que os sı́mbolos “ +”, “∆”, “ o”, “x”, nesta ordem, indicam
os quartis amostrais. Essa imagem confirma a idéia da existência de padrão espacial nos dados, uma vez que existem conglomerados das categorias. Ainda na Figura 2, pode-se observar
o gráfico da coordenada Y pelos dados (canto superior direito) e no canto inferior esquerdo,
gráfico dos dados pela coordenada X e, por último, no canto inferior direito tem-se o histograma da variável.
6
Figura 2: Plot da variável pH da água para análise exploratória
Apresenta-se na Figura 3 o envelope do variograma. Observa-se que existem pontos fora
do mesmo, o que corrobora com as afirmações realizadas a partir da Figura 2, ou seja, existe
dependência espacial.
Figura 3: Gráfico dos envelopes do variograma da variável pH da água
Uma vez que a dependência espacial foi observada, foram modeladas duas funções - exponencial e Matérn - que explicam de forma razoável a variável em estudo. Para ajustar os
modelos por meio do método da máxima verossimilhança, os valores iniciais dos parâmetros
7
foram: σ2 = 0, 9 ∗ Var(dados) = 0, 0833, τ2 = 0, 1 ∗ Var(dados) = 0, 0092 e para φ o valor inicial
foi a máxima distância entre os pontos.
Na Tabela 2, são apresentadas as estimativa de cada um dos parâmetros para os dois modelos, bem como os valores dos critérios de qualidade de ajuste RD, AIC e BIC, para a comparação
dos mesmos.
Tabela 2: Valores dos parâmetros estimados dos modelos ajustado para pH da água
Modelo
σ̂2
τ̂2
φ̂
βˆ0
RD
AIC
BIC
Exponencial 0,0694 0,0392 3,8699 5,0156 36,10 41,197 55,120
Matérn
0,0562 0,0474 2,2218 5,0219 45,75 42,231 56,154
Observando a Tabela 2, nota-se que o modelo melhor ajustado é aquele que utiliza a função
exponencial, pois ao ser comparado com o modelo ajustado com a função Matérn, possui menores valores para AIC e BIC, além de um melhor valor de RD, logo este é o escolhido. Desta
forma, realizou-se a krigagem da área em estudo com o modelo ajustado com a função exponencial, para a variável pH da água (Figura 4).
Figura 4: Krigagem da variável pH da água
Analisando a Figura 4, nota-se que os maiores valores da variável (em verde) estão no canto
da área estudada e ao centro da área os valores são menores (em vermelho).
3.1
Análise bayesiana da Variável pH da água
Para a abordagem bayesiana, ajustou-se apenas o modelo mais bem ajustado por meio da
abordagem frequentista, ou seja, o modelo utilizando a função exponencial. Os resultados da
análise bayesiana dos dados são apresentados na Tabela 3.
8
Tabela 3: Intervalo de Credibilidade para os parâmetros
Parâmetro Mı́nimo 1o Quartil Mediana Média 3o Quartil
βˆ 0
4,779
5,007
5,042
5,039
5,078
σˆ2
0,062
0,083
0,107
0,105
0,124
φ̂
τˆ2
0,821
0
0,821
0
1,641
0
1,289
0
1,641
0
Máximo
5,312
0,342
4,924
0
Observa-se pela Tabela 3 que os valores das estimativas dos parâmetros são muito próximos
dos resultados obtidos pela análise frequentista. A proximidade nas estimativas pode ser atribuı́da
à escolha das funções de distribuição a priori não-informativas para os parâmetros do modelo.
Uma vez ajustado o modelo que utiliza a função exponencial, realizou-se a krigagem (Figura
5)
Figura 5: Krigagem da variável pH da água realizada pela abordagem bayesiana
Observa-se que a krigagem efetuada utilizando-se a abordagem bayesiana (Figura 5), está
muito próxima àquela encontrada do ponto de vista clássico (Figura 4). Como explicitado anteriormente, isso deve-se, principalmente, a escolha de prioris não-informativas para os parâmetros
do modelo em estudo.
4
Conclusões
Por meio do estudo proposto, foi possı́vel avaliar a variabilidade da variável pH da água na
área em estudo, localizada em Unaı́-MG. Para a krigagem foi utilizado o modelo que utiliza
como função de correlação a função exponencial, que se sobressaiu ao modelo com a função
Matérn, por meio dos critérios de seleção RD, AIC e BIC. Além da abordagem frequentista,
9
foi realizada também a abordagem bayesiana, que retornou resultados muito próximos. Essa
proximidade nos resultados possivelmente decorre da utilização de distribuições a priori nãoinformativas.
5
Bibliografia
Referências
[1] CAMBARDELLA, C.A.; MOOMAN, T.B.; NOVAK, J.M.; PARKIN, T.B.; KARLEM,
D.L.; TURVO, R.F.; KONOPA, A.E. Field scale variability of soil properties in central
lowa soil. Soil Sci. Am. J., 47:1501-1511,1994.
[2] DRUCK, S.; CARVALHO,M.S.; CÂMARA,G.; MONTEIRO,A.M.V. Análise espacial
de dados geográficos. Brası́lia, EMBRAPA, 2004 (ISBN: 85-7383-260-6)
[3] KOLMAN, B. Bayesian statistics: an introduction. New York: Hodder Arnold, 2004.
462p.
[4] LANDIM, P.M.B. Sobre Geoestatı́stica e mapas. Terrae Didatica, v. 2, n. 1, p.19-33,
2006. Disponı́vel em: ¡http://www.ige.unicamp.br/terradidatica/¿. Acesso em: 05 nov.
2012.
[5] MATHERON, G. Principles of geostatistics. Economic Geology, Austin, v. 58, p. 12461266, 1963.
[6] R DEVELOPMENT CORE TEAM. R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria, 2011. Disponı́vel
em:
http://www.R-project.org.
[7] RIBEIRO JÚNIOR, P.J.; DIGGLE, P.J. geoR: A package for geoestatistical analysis. RNEWS, Vienna, v.1, n.2, p.14-18, 2001.
[8] RIGHETTO, A. J. Avaliação de modelos geoestatı́sticos multivariados. 2013. 90 p.
Dissertação (Mestrado em Ciências Exatas) - Escola Superior de Agricultura “Luiz de
Queiroz”, Universidade de São Paulo, Piracicaba, 2013.
[9] SALVIANO, A.A.C. Variabilidade de atributos de solo e de Crotalaria juncea em solo
degradado do municı́pio de Piracicaba-SP. 1996. 91 p. Tese (Doutorado em Ciências
Exatas) - Escola Superior de Agricultura “Luiz de Queiroz”, Universidade de São Paulo,
Piracicaba, 1996.
10
[10] VIEIRA, S.R. Geoestatı́stica em estudos de variabilidade espacial do solo. Tópicos em
ciência do solo. Viçosa, v. 1, p. 1-54, 2000.
11