Escola Superior de Agricultura “Luiz de Queiroz” Universidade de São Paulo Dissertação de Mestrado

Bruno H F Fonseca - bhffonse@esalq.usp.br

Estatística e Experimentação Agronômica ESALQ/USP

Orientador: Paulo Justiniano Ribeiro Jr.

Laboratório de Estatística e Geoinformação LEG/UFPR

Piracicaba/2009

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  1 / 70

Um estudo sobre estimação e predição em modelos geoestatísticos bivariados

Motivação

Modelos lineares:

Motivação

Modelos geoestatísticos:

Introdução

Modelos Lineares Univariados Modelos Espaciais

Geoestatística Dados de Área Processos Pontuais

Modelos Geoestatísticos Multivariados

Objetivos

Revisão bibliográfica sobre modelos geoestatísticos bivariados Estudar a estimabilidade paramétrica dos modelos geoestatísticos bivariados Estudar a capacidade de predição espacial dos modelos geoestatísticos bivariados Avaliar a aplicabilidade de modelos geoestatísticos bivariados em dados observacionais

Campos Aleatórios

Uma notação para um campo aleatório Z na localização s_{i} pertecente ao espaço sob estudo G é:

(Z(s_{i}) : s_{i} ∈ G ⊂ R^{d})

Segundo Schmidt e Sansó (2006) e Le e Zidek (2006), a descrição de um campo aleatório é obtida através da distribuição acumulada conjunta finita dimensional de Z(s):

F_{s_{1},s_{2},...,s_{n}}(z_{1}, z_{2}, ..., z_{n}) ≡ P(Z(s_{1}) ≤ z_{1}, Z(s_{2}) ≤ z_{2}, ..., Z(s_{n}) ≤ z_{n})

sendo s = (s_{1}, s_{2}, ..., s_{n}) um conjunto de localizações pertencentes a G.

Campos Aleatórios Gaussianos

Assumindo que o campo aleatório é gaussiano, o processo aleatório é completmente especificado por:

Vetor de média μ = (μ_{1}, μ_{2}, ..., μ_{n})

Matriz de covariâncias

Σ = [ \begin{matrix}σ_{1}^{2} & σ_{1,2} & ⋯ & σ_{1,n} \\ σ_{1,2} & σ_{2}^{2} & ⋯ & σ_{2,n} \\ ⋮ & ⋮ & ... & ⋮ \\ σ_{1,n} & σ_{2,n} & ⋯ & σ_{n}^{2} \end{matrix} ]

Sendo assim a distribuição conjunta de Z(s) é:

f_{s}(z(s)) = (2π)^{−0.5n}|Σ|^{−0.5}exp(−0.5(z(s) − μ)^{t}Σ^{−1}(z(s) − μ))

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  8 / 70

Campos Aleatórios Gaussianos

Em geral não é possível ter mais de uma realização do processo, sendo assim, outras suposições devem ser impostas sobre os campos aleatórios para análises inferenciais:

Estacionariedade Forte

f_{s}(z(s)) = f_{s+h}(z(s + h)), para qualquer incremento h

Isotropia

Dessa forma o campo aleatório é dito ser homogêneo.

Campos Aleatórios Gaussianos

Em alguns casos a suposição de estacionariedade não é válida: Problemas com a média Problemas com a variância

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  10 / 70

Se Z(s) é gaussiano e homogêneo, Σ deve:

ser positiva definida assegurar o comportamento empírico utilizado em geoestatística

Estruturar Σ dessa forma não é trivial, no entanto: Cov(Z(s_{l}), Z(s_{k})) = σ^{2}ρ(Z(s_{l}), Z(s_{k})), para todo (s_{l}, s_{k}) ∈ s existem funções de correlação que asseguram a validade de Σ Funções de correlação válidas:

Família Matérn

ρ(h) = 2^{κ} − Γ(κ)^{−1}(h/ϕ)^{κ}K_{κ}(h/ϕ)

os parâmetros dessa função são ϕ > 0 e κ > 0, que são vinculados a escala com a dimensão de distância e suavidade do processo e K_{κ}(.) é a função Bessel de ordem κ.

Família Exponencial Potência

ρ(h) = exp (h/ϕ)^{κ}

essa família também possui dois parâmetros, com as mesmas interpretações da família Matérn, no entanto agora κ limitado no intervalo [0, 2].

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  12 / 70

Modelo Geoestatístico Univariado

Com as estruturas paramétricas definidas, pode-se pensar na modelagem:

Y (s) = X β + Z(s) + ϵ, s = (s_{1}, s_{2}, ..., s_{n}),  (1)

Y (s) é o vetor dos valores observados da variável de interesse X é a matriz de delineamento β é o vetor de parâmetros associados a média da resposta Z(s) é um campo gaussiano com vetor de médias nulo e Σ válida ϵ é o vetor dos ruídos brancos que são i.i.d. de uma N(0, τ^{2})

Modelo Geoestatístico Univariado

Definida a modelagem para a resposta, tem-se que:

Y (s) ∼ N(μ; Σ_{Y}),

sendo μ = X β e Σ_{Y} = Σ + τ^{2}I , I é uma matriz identidade. Logo θ = (β, σ^{2}, ϕ^{∗}, τ^{2}) é o vetor paramétrico associado a (1). Sendo ϕ^{∗} o vetor de parâmetros associado a função de correlação adotada.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  14 / 70

Estimação em modelo Geoestatístico Univariado

O método da máxima verossimilhança é uma abordagem para estimar θ e consiste em encontrar θ̂ que maximize a função:

l(θ; Y (s)) = −0.5[ln(|Σ_{Y}|) + n ln(2π) + (Y (s) − μ)^{t}Σ_{Y}^{−1}(Y (s) − μ)]

Diggle e Ribeiro Jr (2007) adotam a utilização da reparametrização ν = τ/σ, sendo assim:

θ^{∗} = (β, σ^{2}, ϕ^{∗}, ν^{2})

Σ_{Y} = σ^{2}V

sendo V uma matriz que depende apenas de θ_{c} = (ϕ^{∗}, ν^{2}).

Estimadores de máxima verossimilhança com forma analítica:

β̂ = (X^{t}V^{−1}X )^{−1}(X^{t}V^{−1}Y (s))
σ̂^{2} = n^{−1}(Y (s) − X β)^{t}V^{−1}(Y (s) − X β)

As expressões de β̂ e σ̂^{2} são funções de θ_{c}, além das formas fechadas para os estimadores, é possível encontrar a matriz de informação de Fisher observada para os mesmos e portanto a matriz de covariâncias associada a (β̂; σ̂^{2}) (DIGGLE; RIBEIRO Jr., 2007).

Utilizando β̂ e σ̂^{2} na função de verossimilhança, encontra-se a função de verossimilhança concentrada, que depende apenas de θ_{c}. Os métodos numéricos de estimação podem ser utilizados para encontar:

Estimativa de θ_{c}

Estimativa da matriz Hessiana, denotada por H

Logo, utilizando a propriedade de invariância dos estimadores de máxima verossimilahnça pode-se encontrar β̂, σ̂^{2} e τ̂^{2}.

Utilizando as propriedades assintóticas dos estimadores de máxima verossimilhança e o método delta tem-se que:

θ̂ ∼ N(θ; Σ_{θ})

sendo Σ_{θ} = Δ^{t}Σ_{θ^{∗}}Δ, onde a i-ésima coluna de Δ é o vetor \frac{∂l(θ_{i})}{∂θ^{∗}} e:

Σ_{θ^{∗}} = [ \begin{matrix}Σ_{β,σ^{2}} & O \\ O^{t} & Σ_{θ_{c}} \end{matrix} ]

sendo Σ_{β,σ^{2}} a matriz de covariâncias de (β̂, σ̂^{2}), que possui forma analítica, Σ_{θ_{c}} = −H^{−1} é a matriz de covariâncias de θ̂_{c}, e O uma matriz de zeros, pois não é possível recuperar as covariâncias entre (β̂, σ̂^{2}) e θ̂_{c}.

Krigagem em modelo Geoestatístico Univariado

Processo de predição do campo aleatório em localizações não amostradas, supondo que o interesse seja fazer predição de Z em s^{∗}, tem-se que a distribuição de (Z(s^{∗})|Y (s)) gera as predições com melhor precisão:

E(Z(s^{∗})|Y (s)) = μ_{Z} + Σ_{Z,Y}Σ_{Y}^{−1}(Y (s) − X β)
Var(Z(s^{∗})|Y (s)) = Σ_{Z} − Σ_{Z,Y}Σ_{Y}^{−1}Σ_{Z,Y}^{t}

sendo μ_{Z} a média de Z(s^{∗}), Σ_{Z,Y} é a matriz de covariâncias cruzadas entre Z(s^{∗}) e Y (s) e Σ_{Z} a matriz de covariâncias de Z(s^{∗}).

Modelos Geoestatísticos Biivariados

Quando existem dois campos aleatórios de interesse, pode-se pensar na modelagem:

Y_{i} = μ_{i} + Z_{i}, i = 1, 2  (2)

Y_{1} e Y_{2} são os vetores dos valores observados das respostas μ_{i} é o vetor de médias associado a Y_{i} e pode incluir covariáveis Z_{i} é um campo gaussiano com vetor de médias nulo e Σ_{i} válida A distribuição de Y = (Y_{1}, Y_{2}) é de interesse, Y ∼ N(μ, Σ_{Y}), sendo μ = (μ_{1}, μ_{2}) e:

Σ_{Y} = [ \begin{matrix}Σ_{1} & Σ_{1,2} \\ Σ_{1,2}^{t} & Σ_{2} \end{matrix} ]

sendo Σ_{i} a matriz de covariância de Y_{i}, i = 1, 2 e Σ_{1,2} a matriz de covariâncias cruzadas entre as respostas. A especificação do modelo deve assegurar que Σ_{Y} seja válida:

Modelo gaussiano bivariado com componente de correlação parcialmente comum - BGCCM Modelo bivariado de co-regionalização - BCRM A distribuição de Y = (Y_{1}, Y_{2}) é de interesse, Y ∼ N(μ, Σ_{Y}), sendo μ = (μ_{1}, μ_{2}) e:

Σ_{Y} = [ \begin{matrix}Σ_{1} & Σ_{1,2} \\ Σ_{1,2}^{t} & Σ_{2} \end{matrix} ]

sendo Σ_{i} a matriz de covariância de Y_{i}, i = 1, 2 e Σ_{1,2} a matriz de covariâncias cruzadas entre as respostas. A especificação do modelo deve assegurar que Σ_{Y} seja válida:

Modelo gaussiano bivariado com componente de correlação parcialmente comum Modelo bivariado de co-regionalização

BGCCM

Diggle e Ribeiro Jr. (2007) propuseram um modelo que decompõe os termos latentes de (2):

Z_{i} = σ_{0i}S_{0} + σiS_{i}, i = 1, 2,  (3)

σ^{∗} = (σ_{01}, σ_{1}, σ_{02}, σ_{2}) um vetor de parâmetros de dispersão S_{0}, S_{1} e S_{2} são campos aleatórios gaussianos mutuamente independentes, com vetores de médias nulos, variâncias unitárias e correlações determinadas por funções de correlação válidas

BGCCM

Com as decomposições de (3) tem-se que (2) fica na forma:

\{ \begin{matrix}Y_{1} = μ_{1} + σ_{01}S_{0} + σ_{1}S_{1} \\ Y_{2} = μ_{2} + σ_{02}S_{0} + σ_{2}S_{2} \end{matrix}

Utilizando propriedades básicas de covariância tem-se:

Var(Y_{1}) = σ_{01}^{2} + σ_{1}^{2} Var(Y_{2}) = σ_{02}^{2} + σ_{2}^{2} Cov(Y_{1,l}, Y_{1,k}) = σ_{01}^{2}ρ_{0}(h) + σ_{1}^{2}ρ_{1}(h) Cov(Y_{2,l}, Y_{2,k}) = σ_{02}^{2}ρ_{0}(h) + σ_{2}^{2}ρ_{2}(h) Cov(Y_{1,l}, Y_{2,k}) = σ_{01}σ_{02}ρ_{0}(h)

sendo que os índices l e k referessem à duas localizações espaciais quaisquer sepradas pela distância h.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  24 / 70

BGCCM

Sendo assim, a distribuição de probabilidade de Y = (Y_{1}, Y_{2}) está completamente estruturada e depende de:

θ = (β^{∗}, σ^{∗}, ϕ_{0}^{∗}, ϕ_{1}^{∗}, ϕ_{2}^{∗})

Os processos de estimação e predição seguem a mesma tecnicalidade aplicada nos modelos univariados, porém existe a necessidade de um número maior de reparametrizações.

BCRM

Gelfand et al. (2004) propôs outra abordagem para decompor os termos latentes de (2):

\{ \begin{matrix}Z_{1} = σ_{11}S_{1} \\ Z_{2} = σ_{12}S_{1} + σ_{22}S_{2} \end{matrix}

σ^{∗} = (σ_{11}, σ_{12}, σ_{22}) um vetor de parâmetros de dispersão S_{1} e S_{2} são campos aleatórios gaussianos mutuamente independentes, com vetores de médias nulos, variâncias unitárias e correlações determinadas por funções de correlação válidas

BCRM

Com as decomposições utilizadas no BCRM tem-se:

\{ \begin{matrix}Y_{1} = μ_{1} + σ_{11}S_{1} \\ Y_{2} = μ_{2} + σ_{12}S_{1} + σ_{22}S_{2} \end{matrix}

Utilizando propriedades básicas de covariância tem-se:

Var(Y_{1}) = σ_{11}^{2} Var(Y_{2}) = σ_{12}^{2} + σ_{22}^{2} Cov(Y_{1,l}, Y_{1,k}) = σ_{11}^{2}ρ_{1}(h) Cov(Y_{2,l}, Y_{2,k}) = σ_{12}^{2}ρ_{1}(h) + σ_{22}^{2}ρ_{2}(h) Cov(Y_{1,l}, Y_{2,k}) = σ_{11}σ_{12}ρ_{1}(h)

sendo que os índices l e k referessem à duas localizações espaciais quaisquer sepradas pela distância h.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  27 / 70

BCRM

Sendo assim, a distribuição de probabilidade de Y = (Y_{1}, Y_{2}) está completamente estruturada e depende de:

θ = (β^{∗}, σ^{∗}, ϕ_{1}^{∗}, ϕ_{2}^{∗})

Os processos de estimação e predição seguem a mesma tecnicalidade aplicada nos modelos univariados, porém existe a necessidade de um número maior de reparametrizações.

Material e Métodos

Com intuito de comparar as capacidades de estimação sob o enfoque frequentista e de predição dos modelos bivariados, a pesquisa foi dividida em duas etapas:

Estudo de Simulação Análise de Dados de Solos

Cabe ressaltar que todas análises foram conduzidas no ambiente R de programação (R Development Core Team, 2006), utilizando o pacote geoR (Ribeiro Jr. e Diggle, 2001). No entanto, para os modelos bivariados a maioria dos métodos utilizados não está implementada nos pacotes do R.

Material e Métodos - Simulação

No estudo de simulação foi lançado mão das seguintes características nas duas abordagens de modelos bivariados:

Três configurações amostrais de localizações espaciais

n_{1} = n_{2} = 100 e dados completamente co-locados n_{1} = 100 e n_{2} = 50 e dados completamente co-locados n_{1}=n_{2} = 100 e apenas metade dos dados co-locados

Simulação de 1000 conjuntos de dados em cada configuração Divisão dos dados simulados em um grupo de controle e um grupo de modelagem

Material e Métodos - Simulação

Na simulação do BGCCM foram utilizadas as seguintes configurações paramétricas:

μ_{1} = 150 e μ_{2} = 60 σ_{01} = 8, σ_{1} = 4, σ_{02} = 5 e σ_{2} = 2 Função de Matèrn foi adotada para as três funções de correlação

ϕ_{0} = 0, 25, ϕ_{1} = 0, 2 e ϕ_{2} = 0, 2

κ_{0} = 0, 5, κ_{1} = 0, 5 e κ_{2} = 0, 5
Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  31 / 70

Material e Métodos - Simulação

Na simulação do BCRM foram levadas em consideração características de variabilidade e correlação parecidas com as do BGCCM:

μ_{1} = 150 e μ_{2} = 60 σ_{11} = 9, σ_{12} = 5 e σ_{22} = 2 Função de Matèrn foi adotada para as três funções de correlação

ϕ_{1} = 0, 25 e ϕ_{2} = 0, 2

κ_{1} = 0, 5 e κ_{2} = 0, 5
Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  32 / 70

Para calcular as estimativas por máxima verossimilhança utilizou-se as seguintes reparametrizações:

BGCCM

σ = σ_{01}

η = σ_{02}/σ_{01}
ν_{1} = σ_{1}/σ_{01}
ν_{2} = σ_{2}/σ_{01}

BCRM

σ = σ_{11}

ν_{1} = σ_{12}/σ_{11}
ν_{2} = σ_{22}/σ_{11}

Com essas reparametrizações é possível encontrar a função de máxima verossimilhança concentrada. Logo existe forma analítica para os parâmetros de média e para σ, os demais parâmetros foram estimados pelo método de Nelder e Mead (1965). E novamente as propriedades assintóticas e de invariância foram utilizadas para encontrar a distribuição de probabilidade do vetor de parâmetrso.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  33 / 70

Para avaliar a qualidade das estimativas foram analisados:

erros absolutos médios erros quadráticos médios

Cada vetor de parâmetros estimados foram utilizados para fazer a predição nas localizações do grupo de controle:

intervalos de confiança marginais de predição erros absolutos médios erros quadráticos médios

Material e Métodos - Dados de Solos

Pesquisa realizada em 2006 Propriedade agrícola com 51,8 ha Solo argissolo com textura média Dois históricos de manejo 67 localizações amostradas com grid regular Coordenadas no sistema UTM Variáveis químicas analisadas: pH e saturação por bases

Resultados e Discussões - Simulação

θ mínimo 1^{o} quartil mediana média 3^{o} quartil máximo
β_{01} 137, 03 147, 39 150, 01 149, 99 152, 48 161, 82
β_{02} 52, 18  58, 56  60, 02 59, 98  61, 43  66, 32
σ_{01} 4, 13  6, 24  7, 22  7, 41  8, 35  16, 47
σ_{1} > 0, 001  2, 14  3, 58  3, 18  4, 29  7, 81
σ_{02} 2, 12  3, 82  4, 44  4, 58  5, 17  9, 11
σ_{2} > 0.001  0, 38  1, 65  1, 54  2, 40  5, 20
ϕ_{0}  0, 05  0, 14  0, 19  0, 22  0, 27  1, 00
ϕ_{1} > 0.001  0, 11  0, 15  0, 19  0, 22  1, 00
ϕ_{2} > 0.001  0, 11  0, 17  0, 21  0, 28  1, 00

Tabela: Estatísticas descritivas das estimativas por máxima verossimilhança dos paramêtros do BGCCM, usando as amostras simuladas do BGCCM com dados co-locados e balanceados

θ mínimo 1^{o} quartil mediana média 3^{o} quartil máximo
β_{01} 0, 78  4, 97  8, 58 12, 03 15, 10  137, 48
β_{02} 0, 40  1, 81  3, 06  4, 40  5, 55  43, 38
σ_{01} 0, 04  0, 10  0, 13  0, 14  0, 17  0, 68
σ_{1} > 0, 001  0, 04  1, 93 28, 97  9, 63  2467, 42
σ_{02} 0, 01  0, 07  0, 33  4, 78  1, 67  386, 51
σ_{2} > 0.001  0, 12  0, 42  6, 93  1, 92  797, 07
ϕ_{0}  0, 00  0, 00  0, 002 0, 01  0, 009  2, 48
ϕ_{1}  0, 00  0, 00  0, 003 0, 04  0, 01  2, 56
ϕ_{2}  0, 00  0, 00  0, 002 0, 07  0, 02  7, 70

Tabela: Estatísticas descritivas das variâncias das estimativas por máxima verossimilhança dos paramêtros do BGCCM, usando as amostras simuladas do BGCCM com dados co-locados e balanceados

Erro β_{01} β_{02} σ_{01} σ_{1} σ_{02} σ_{2}  ϕ_{0}  ϕ_{1}  ϕ_{2}
AM 2, 98 1, 71 1, 44 1, 41 0, 93 1, 01 0, 094 0, 102 0, 116
QM 14, 04 4, 79 3, 17 3, 65 1, 27 1, 52 0, 015 0, 021 0, 024

Tabela: Erros absolutos médios e erros quadráticos médios das estimativas por máxima verossimilhança dos paramêtros do BGCCM, usando as amostras simuladas do BGCCM com dados co-locados e balanceados

θ  mínimo 1^{o} quartil mediana média 3^{o} quartil máximo
EAM_{Y_{1}} 1, 79  3, 29  3, 78  3, 83  4, 30  6, 21
EAM_{Y_{2}} 1, 14  2, 00  2, 30  2, 31  2, 58  4, 06
EQM_{Y_{1}} 5, 67  17, 29  22, 24 23, 56 28, 55  54, 45
EQM_{Y_{2}} 2, 05  6, 32  8, 13  8, 50  10, 20  21, 95

Tabela: Estatísticas descritivas dos erros de krigagens com as estimativas paramétricas do BGCCM, usando as amostras simuladas do BGCCM com dados co-locados e balanceados

θ mínimo 1^{o} quartil mediana média 3^{o} quartil máximo
\overline{EP}_{Y_{1}} 3, 70  4, 52  4, 74  4, 75  4, 99  5, 83
\overline{EP}_{Y_{2}} 2, 21  2, 71  2, 84  2, 85  2, 98  3, 68

Tabela: Estatísticas descritivas das médias dos erros padrões de krigagens com as estimativas paramétricas do BGCCM, usando as amostras simuladas do BGCCM com dados co-locados e balanceados

PARA NÃO CAIR NA REPETIÇÃO DE TABELAS E RESULTADOS, COMENTAREI AQUI DE UMA FORMA MAIS GERAL OS RESULTADOS DE TODO O ESTUDO DE SIMULAÇÃO

Resultados e Discussões - Dados de Solos

Figura: Gráficos de círculos - o gráfico a esquerda é para a variável pH e o da direita é da saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  42 / 70

Resultados e Discussões - Dados de Solos

Variável Mediana Média D.P.
Saturação 56, 00 53, 27 10, 05
pH  5, 00  4, 94 0, 42

Tabela: Médias, medianas e desvios padrões do pH e da saturação por bases Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  43 / 70

Resultados e Discussões - Dados de Solos

Figura: Gráficos descritivos do padrão espacial do pH

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  44 / 70

Resultados e Discussões - Dados de Solos

Figura: Gráficos descritivos do padrão espacial da saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  45 / 70

Resultados e Discussões - Dados de Solos

Figura: Gráficos de possíveis transformações de variáveis de Box-Cox, o gráfico da esquerda é relativo ao pH e o gráfico da direita é relativo à saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  46 / 70

\begin{matrix}β & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ 4, 907 & 0, 077 & 0, 163 & 560 & 0, 5 & −25, 49 & 58, 98 \\ 4, 905 & 0, 101 & 0, 149 & 510, 6 & 1 & −25, 69 & 59, 38 \\ 4, 904 & 0, 106 & 0, 251 & 650, 5 & 1, 5 & −25, 57 & 59, 15 \\ 4, 903 & 0, 106 & 0, 288 & 569, 7 & 2 & −25, 45 & 58, 90 \\ 4, 903 & 0, 106 & 0, 297 & 491, 1 & 2, 5 & −25, 36 & 58, 73 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para o pH com média constantes

\begin{matrix}β_{0} & β_{1} & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ 4, 725 & 0, 397 & 0 & 0, 125 & 49, 19 & 0, 5 & −23, 44 & 56, 89 \\ 4, 723 & 0, 402 & 0 & 0, 125 & 33, 42 & 1 & −23, 35 & 56, 70 \\ 4, 723 & 0, 403 & 0, 014 & 0, 111 & 28, 43 & 1, 5 & −23, 33 & 56, 66 \\ 4, 779 & 0, 227 & 0, 110 & 0, 061 & 349, 3 & 2 & −24, 86 & 59, 73 \\ 4, 786 & 0, 204 & 0, 109 & 0, 086 & 353, 1 & 2, 5 & −24, 83 & 59, 65 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para o pH com média induzida pela área de manejo

Resultados e Discussões - Dados de Solos

\begin{matrix}β_{0} & β_{1} & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ −607, 48 & 0, 001 & 0 & 0, 111 & 49, 32 & 0, 5 & −19, 66 & 49, 31 \\ −608, 51 & 0, 001 & 0 & 0, 112 & 34, 28 & 1 & −19, 44 & 48, 89 \\ −608, 33 & 0, 001 & 0 & 0, 112 & 27, 58 & 1, 5 & −19, 39 & 48, 77 \\ −608, 36 & 0, 001 & 0, 008 & 0, 104 & 24, 53 & 2 & −19, 36 & 48, 73 \\ −608, 46 & 0, 001 & 0, 015 & 0, 097 & 22, 46 & 2, 5 & −19, 35 & 48, 70 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para o pH com média induzida pelas coordenadas oeste-leste

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  49 / 70

\begin{matrix}β & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ 49, 53 & 47, 66 & 105, 33 & 700, 00 & 0, 5 & −239, 53 & 487, 05 \\ 48, 53 & 59, 57 & 120, 62 & 625, 58 & 1 & −239, 76 & 487, 51 \\ 47, 98 & 62, 40 & 124, 31 & 516, 78 & 1, 5 & −239, 79 & 487, 57 \\ 47, 27 & 63, 29 & 149, 69 & 489, 58 & 2 & −239, 75 & 487, 50 \\ 47, 18 & 63, 51 & 139, 18 & 400, 00 & 2, 5 & −239, 70 & 487, 40 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para a saturação por bases com média constantes

\begin{matrix}β_{0} & β_{1} & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ 47, 91 & 8, 80 & 24, 35 & 53, 68 & 73, 68 & 0, 5 & −238, 05 & 486, 10 \\ 47, 92 & 8, 83 & 38, 78 & 39, 28 & 59, 18 & 1 & −238, 01 & 486, 01 \\ 47, 92 & 8, 86 & 43, 35 & 34, 69 & 50, 33 & 1, 5 & −237, 98 & 485, 96 \\ 47, 92 & 8, 88 & 45, 53 & 32, 47 & 44, 34 & 2 & −237, 96 & 485, 92 \\ 47, 93 & 8, 89 & 46, 79 & 31, 18 & 39, 98 & 2, 5 & −237, 95 & 485, 90 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para a saturação por bases com média induzida pela área de manejo

\begin{matrix}β_{0} & β_{1} & τ^{2} & σ^{2} & ϕ & κ & l(θ) & AIC \\ −14297, 01 & 0, 025 & 8, 98 & 59, 98 & 54, 51 & 0, 5 & −234, 75 & 479, 50 \\ −14300, 66 & 0, 025 & 27, 56 & 41, 46 & 45, 42 & 1 & −234, 72 & 479, 44 \\ −14300, 92 & 0, 025 & 33, 53 & 35, 53 & 39, 51 & 1, 5 & −234, 70 & 479, 39 \\ −14300, 17 & 0, 025 & 36, 39 & 32, 68 & 35, 33 & 2 & −234, 68 & 479, 36 \\ −14298, 98 & 0, 025 & 38, 06 & 31, 03 & 32, 19 & 2, 5 & −234, 67 & 479, 34 \end{matrix}

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao modelo para a saturação por bases com média induzida pelas coordenadas oeste-leste

Resultados e Discussões - Dados de Solos

Utilizando o AIC e os máximos estimados das funções de vereossimilhança tem-se os modelos univariados finais:

Média induzida pelas coordenadas oeste-leste

parâmetros de suavidade da função de correlação igual a 2,5

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  53 / 70

Resultados e Discussões - Dados de Solos

Figura: Gráficos de predições espacias, o gráfico da esquerda é relativo ao pH e o gráfico da direita é relativo à saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  54 / 70

Resultados e Discussões - Dados de Solos

Figura: Densidade e gráfico de quartis dos resdíduos para o modelo final do pH

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  55 / 70

Resultados e Discussões - Dados de Solos

Figura: Densidade e gráfico de quartis dos resdíduos para o modelo final da saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  56 / 70

Resultados e Discussões - Dados de Solos

θ Constantes Oeste-leste Área de manejo
β_{01}  52, 57  1, 54e − 6  47, 83
β_{02}  4, 92  1, 45e − 7  4, 72
β_{1}  −  8, 98e − 5  9, 46
β_{2}  −  8, 41e − 6  0, 40
σ_{01}  8, 51  8, 35  6, 85
σ_{1}  3, 47  3, 46  3, 44
σ_{02}  0, 38  0, 38  0, 31
σ_{2}  7e − 5  5e − 4  1e − 5
ϕ_{0}  59, 37  57, 79  26, 73
ϕ_{1}  46, 35  48, 15  49, 08
ϕ_{2}  66, 82  71, 59  86, 95
κ_{0}  1  1  1, 5
κ_{1}  0, 5  0, 5  0, 5
κ_{2}  0, 5  0, 5  0, 5

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao BGCCM para a saturação por bases e o pH

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  57 / 70

θ Constantes Oeste-leste Área de Manejo
β_{01}  51, 76  1, 51e − 6  47, 93
β_{02}  4, 89  1, 43e − 7  4, 72
β_{1}  −  8, 85e − 5  9, 12
β_{2}  −  8, 36e − 6  0, 40
σ_{11}  8, 26  99, 53  7, 82
σ_{12}  0, 31  3, 67  0, 28
σ_{22}  0, 13  1, 54  0, 14
ϕ_{1}  100, 63  100, 21  53, 32
ϕ_{2}  21, 43  21, 38  19, 34
κ_{1}  0, 5  0, 5  0, 5
κ_{2}  2, 5  2, 5  2, 5

Tabela: Estimativas de máxima verossimilhança dos parâmetros associados ao BCRM para a saturação por bases e o pH

Resultados e Discussões - Dados de Solos

Figura: Gráficos de predições espacias com o BGCCM, o gráfico da esquerda é relativo ao pH e o gráfico da direita é relativo à saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  59 / 70

Resultados e Discussões - Dados de Solos

Figura: Gráficos de predições espacias com o BCRM, o gráfico da esquerda é relativo ao pH e o gráfico da direita é relativo à saturação por bases

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  60 / 70

Obs. Univ. BGCCM BCRM Obs. Univ. BGCCM BCRM
31 37, 98 36, 70 36, 41 47 46, 37 47, 39 46, 98
52 50, 32 48, 47 48, 75 67 58, 89 70, 08 69, 86
41 38, 73 43, 29 42, 70 57 54, 88 59, 06 58, 96
51 42, 33 50, 60 50, 67 56 57, 62 52, 66 52, 93
41 56, 76 39, 96 40, 90 59 60, 97 54, 68 55, 20
43 53, 04 46, 20 46, 51 57 56, 83 56, 93 57, 17
39 44, 37 43, 18 42, 65 52 57, 79 52, 64 52, 95
51 56, 27 48, 48 48, 80 61 56, 38 61, 32 61, 35
64 46, 88 68, 33 68, 03 59 56, 07 54, 46 54, 30
53 47, 06 53, 68 52, 72 54 62, 32 56, 91 57, 16

Tabela: Valores observados da saturação por bases e não utilizados na modelagem e suas respectivas predições utilizando os modelo univariado, o BGCCM e o BCRM

Estatística  Univariado BGCCM BCRM
Média  −0, 34  −0, 50 −0, 50
Desvio Padrão  7, 52  3, 02  2, 84

Tabela: Médias e desvios padrões dos erros de krigagem da saturação por bases nas 20 localizações omitidas no processo de estimação

Considerações Finais

FORMALIZAR AINDA

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  63 / 70

Referências

AZZALINI, A. Statistical inference: based on the likelihood. London: Chapman and Hall, 1996. 341 p.

BANERJEE, S.; CARLIN B.P.; GELFAND, A.E. Hierarchical modelling and analysis for spatial data. Boca Raton: Chapman and Hall, 2004. 453 p.

BANERJEE, S.; GELFAND, A.E. Predict, Interpolation and regression for spatial misaligned data points. Sankhya, Kolkata, v.64, p.227-245, 2002.

BICKEL, P.J.; DOKSUM K.A. Mathematical statistics: basic ideas and selected topics. San Francisco: Holden-Day, 1977. 492 p.

BOGNOLA, I. A.; RIBEIRO JR, P. J.; SILVA, E. A. A.; LINGNAU, C.; HIGA, A. R. Modelagem uni e bivariada da variabilidade espacial de rendimento de Pinus taeda L.. Floresta, Curitiba, v. 38, p. 373-385, 2008.

BROWN, P.J.; LE, N.D.; ZIDEK, J.V. Multivariate spatial interpolation and exposure to air pollutants. The Canadian Journal of Statistics, Vancouver, v.22, p.489-509, 1994.

CHRISTENSEN, W.F.; AMEMIYA, Y. Latent variable analysis of multivariate spatial data. Journal of the American Statistical Association, Alexandria, v.97, p.302-317, 2002.

CHRISTENSEN, W.F.; DIGGLE, P.J.; RIBEIRO Jr., P.J. Analysing positive-valued spatial data: the transformed gaussian model. Geostatistics for environmental applications, Boston, v.11, p.287-298, 2001.

DIGGLE, P.J.; RIBEIRO Jr., P.J. Bayesian inference in Gaussian model based geostatistics. Geographical and Environmental Modelling, New Jersey, v.6, p.129-146, 2002. DIGGLE, P.J.; RIBEIRO Jr., P.J. Model-Based geostatistics. New York: Springer, 2007. 230p.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  65 / 70

FARACO, M.A.; URIBE-OPAZO, M.A.; SILVA E.A.A.; JOHANN J.A.; BORSSOI J.A. Nonstationary multivariate process modeling through spatially varying coregionalization. Revista Brasileira de Ciências do Solo, Viçosa, v.32, p.463-476, 2008.

FERNANDES, M.V.M. Modelos para Processos Espaço-Temporais Inflacionados de Zeros, 2006. 128 p. Dissertação (Mestrado em Estatística) - Instituto de Matemática da Univeridade Federal do Rio de Janeiro, 2006.

FINLEY, A.O.; BANERJEE S.; CARLIN B.P. spBayes: an R package for univariate and multivariate hierarchical point-referenced spatial models. Journal of Statistical Software, Alexandria, v. 19, p.1-24, 2007.

FUENTES, M.; SMITH, R.L. A new class of stationary spatial models. North Caroline: Department of Statistics, North Caroline State University, 2001. Technical (Report, 2534).

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  66 / 70

GELFAND, A.E.; SCHMIDT, A.M.; BANERJEE S.; SIRMANS, C.F. Nonstationary multivariate process modeling through spatially varying coregionalization. Test, Valencia, v.13, p.263-312, 2004.

HANDCOCK, M.; STEIN, M. A Bayesian analysis of kriging. Technometrics, Alexandria, v.35, p.403-410, 1993.

HIGDON, D. Quantitative methods for current environmental issues. Chichester: Wiley, 2002. 185 p.

KITANIDIS, P.K. Parameter uncertainty in estimation of spatial functions: Bayesian analysis. Water Resources Research, New York, v.22 p.499-507, 1986.

LE, D.N.; ZIDEK, J.V. Statistical analysis of environmental space-time processes. New York: Springer, 2006. 327p.

MATÈRN, B. Spatial variation. Verlag, Berlin: Springer, 1986. 365 p.

Bruno H F Fonseca (ESALQ/USP)  11 de Fevereiro de 2009  67 / 70

NELDER, J.A.; MEAD R. A simplex method for function minimization. The Computer Journal, Oxford, v.7: p.308-313, 1965.

RIBEIRO Jr, P.J.; DIGGLE, P.J. geoR: A package for geostatistical analysis. R-NEWS, Vienna, v.1, n.2, p.14-18, 2001. ISSN 1609-3631.

ROYLE, A.J.; BERLINER, M.L. A hierarchical approach to multivariate spatial modeling and prediction. Journal of Agricultural, Biological, and Environmental Statistics, Alexandria, v.4, p.29-56, 1999.

SAMPSON, P.D.; GUTTORP, P. Nonparametric estimation of nonstationary spatial covariance structure. Journal of American Statistical Association, Alexandria, v.87, p.108-119, 1992.

SCHABENBERGER, O.; GOTWAY, C.A. Statistical methods for spatial data analysis. Boca Raton: Chapman and Hall / CRC, 2005. 488p.

SCHMIDT, A.M.; GELFAND, A.E. A Bayesian coregionalization approach for multivariate pollutant data. Journal of Geophysical Research - Atmospheres, Florida, v.108, p.8783, 2003.

SCHMIDT, A.M.; SANSÓ, B. Modelagem bayesiana da estrutura de covariância de processos espaciais e espaço-temporais. In: SIMPÓSIO NACIONAL DE PROBABILIDADE E ESTATÍSTICA, 14, 2006. Caxambú, Minicurso... São Paulo: Associação Brasileira de Estatística, 2006. 151 p.

SCHMIDT, A.M.; O’HAGAN, A. Bayesian inference for nonstationary spatial covariance structure via spatial deformations, Journal of Royal Statistical Society, Oxford, v.65, p.743-758, 2003.

SCHULLER, J.K. A review and integrating analysis of spatially-variable crop control of crop production. Fertilizer Research, Amsterdam, v.33, p.1-34, 1992.

SILVA, A.S. Modelos gaussianos geoestatísticos espaço-temporais e aplicações, 2006. 70 p. Dissertação (Mestrado em Estatística e Experimentação Agronômica) - Escola Superior de Agricultura ”Luiz de Queiroz” da Univeridade de São Paulo, 2006.

SUN, W.; LE, N.D.; ZIDEK, J.V.; BURNETT, R. Assessment of a Bayesian multivariate interpolation approach for health impact studies. Environmetrics, Washington, v.9, p.565-586, 1998.

R: A Language and Environmental for Statistical Computing, R Development Core Team, R Foundation for Statistical Computing, Vienna, Austria, 2006, ISBN 3-900051-07-0. Disponível em: ¡http://www.R-project.org¿. Acesso em 28 maio 2007.