{\setlength{\baselineskip}{1.0\baselineskip}

\section{Material e Métodos}
Alguns trabalhos recentemente publicados vêm utilizando modelos geoestatísticos bivariados para estruturar a dependência espacial entre e dentre campos aleatórios de interesse (SCHMIDT E SANSÓ 2006; BOGNOLA ET AL. 2008). 
No entanto, não existem pesquisas comparando o BGCCM e o BCRM sob o enfoque frequentista de estimação. Esse trabalho concentra-se em utilizar tais modelos e verificar suas 
características, vantagens e desvantagens na estimação de parâmetros e predição 
espacial em análise de dados geoestatísticos.

Com o intuito de verificar a qualidade das estimativas por máxima verossimilhança e das krigagens utilizando os modelos geoestatísticos bivariados, inicialmente foi conduzido um estudo de simulação com diversas configurações para o BGCCM e o BCRM.

Após o estudo de simulação foram analisados dados observacionais provinientes de uma pesquisa sobre qualidade do solo de uma propriedade agrícola.

Todas as análises e resultados foram obtidos através do 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, sendo assim, consta no apêndice as programações utilizadas.

\subsection{Estudo de Simulação}
Essa fase do trabalho é muito importante para detectar possíveis dificuldades no ajuste de modelos geoestatísticos bivariados e verificar as vantagens e desvantagens de cada abordagem considerada. Nesse contexto, foram selecionadas 3 configurações amostrais de localizações espaciais para cada abordagem e foram simulados 1000 conjuntos de dados dentro de um quadrado de lado igual a 1 e com uma malha irregular de localizações para cada configuração. Para avaliar as precisões das estimativas e das predições espaciais, cada vetor $Y=(Y_{1},Y_{2})$ de dados simulados foi dividido em um grupo de modelagem e um grupo de controle. O grupo de modelagem foi utilizado para ajustar os modelos bivariados e o grupo de controle, que possui 20 observações de cada resposta e suas respectivas localizações espaciais, foi utilizado para avaliar os erros de predição nas localizações desse grupo.

Para simular observações com o BGCCM, foram considerados os parâmetros de média e variabilidade sempre fixos em $\mu_{1}=150$, $\mu_{2}=60$, $\sigma_{01}=8$, $\sigma_{1}=4$, $\sigma_{02}=5$ e $\sigma_{2}=2$, logo, $Var(Y_{1})=\sigma_{01}^2 + \sigma_{1}^2 = 80$, $Var(Y_{1})=\sigma_{02}^2 + \sigma_{2}^2=29$. A função de Matèrn foi adotada para as três funções de correlação com valores de parâmetros $\phi_{0}=0,25$, $\phi_{1}=0,2$, $\phi_{2}=0,2$, $\kappa_{0}=0,5$, $\kappa_{1}=0,5$ e $\kappa_{2}=0,5$. Utilizando essa configuração paramétrica tem-se que os dois campos aleatórios são homogêneos e os parâmetros escolhidos para estruturar as correlações geram alcances práticos, distância $h$ em que as $\rho(h)=0.05$, aproximadamente iguais a $0,75$, $0,60$ e $0,60$ para $\rho_{0}$, $\rho_{1}$ e $\rho_{2}$, respectivamente. Já para a amostragem de localizações, foram consideradas três configurações distintas: na primeira foram utilizados $n_{1}=n_{2}=100$ e as localizações das duas respostas completamente coincidentes, dados co-locados; na segunda foram utilizados dados completamente co-locados, porém desbalanceados, sendo $n_{1}=100$ e $n_{2}=50$; e na última configuração foi considerado $n_{1}$=$n_{2}=100$ porém apenas metade dos dados co-locados.

Com o BCRM, foram consideradas as configurações das localizações e dos parâmetros de média idênticas as do BGCCM. Os demais parâmetros foram fixados em $\sigma_{1}=9$, $\sigma_{12}=5$, $\sigma_{2}=2$, $\phi_{1}=0,25$, $\phi_{2}=0,2$, $\kappa_{1}=0,5$ e $\kappa_{2}=0,5$. Tais configurações foram utilizadas para gerar características similares em ambos os casos permitindo
comparações entre os modelos considerados.

Os parâmetros foram estimados por máxima verossimilhança em cada uma das observações de $Y$ utilizando o mesmo modelo do qual a observação é proviniente. Para facilitar a estimação de $\theta$ foram utilizadas reparametrizações. No BGCCM foram definidos $\sigma=\sigma_{01}$,  $\eta=\sigma_{02}/\sigma_{01}$, $\nu_{1}=\sigma_{1}/\sigma_{01}$ e $\nu_{2}=\sigma_{2}/\sigma_{01}$, já no BCRM foram utilizados $\sigma=\sigma_{11}$, $\nu_{1}=\sigma_{12}/\sigma_{11}$ e $\nu_{2}=\sigma_{22}/\sigma_{11}$. Nas duas abordagens pode-se escrever $\Sigma_{Y}=\sigma^{2}V$, sendo assim, existe forma analítica para $\widehat{\mu}$ e $\widehat{\sigma}$ e $V$ é função dos demais parâmetros, os quais foram estimados maximizando a função de máxima verossimilhança concentrada pelo método de Nelder e Mead (1965).

As tecnicalidades para encontrar $\widehat{\theta}$ e sua distribuição de probabilidade assintótica são similares as utilizadas nos modelos geoestatísticos univariados, onde foi utilizada a propriedade de invariância para calcular todas as estimativas e o método delta para encontrar a matriz de variância de $\widehat{\theta}$. Cabe ressaltar que em nenhum caso foram estimados os parâmetros de suavidade das funções de correlação que foram considerados fixos.

Para avaliar a qualidade das estimativas, foram utilizados erros absolutos médios e erros quadráticos médios.

Cada vetor de parâmetros estimados foi utilizado para fazer as krigagens dos campos aleatórios nas localizações do grupo de controle. Sendo assim, para cada modelo ajustado, existem 20 valores observados de cada resposta e suas predições espaciais. Para avaliar a precisão dessas predições, foram calculados intervalos de confiança marginais para as predições, erros quadráticos médios e erros absolutos médios.

\subsection{Dados sobre a qualidade do solo}
Os dados analisados com as técnicas apresentadas foram obtidos por meio de uma pesquisa realizada em julho 2006 na fazenda Tupã, localizada no município de Echaporan/SP, que possui 51,8ha de área, solo argissolo de textura média e dois históricos de manejo distintos, soja numa região e pastagem na outra. Para analisar as variáveis químicas do solo, foram amostradas 67 localizações, com grid regular a cada hectare, georeferenciadas no sistema Universal Transverse Mercator (UTM). Com a utilização de um aparelho GPS foram medidos os atributos de interesse.

Dois parâmetros químicos do solo levantados pelos pesquisadores serão
modelados nesse trabalho: saturação por bases e pH do solo, a primeira
variável é uma medida de capacidade do solo reter bons nutrientes N,
P, K, Ca, Mg e a segunda variável mede a acidez do solo, sendo que, existe uma forte correlação entre essas duas respostas, o
que justifica a tentativa de modelagem bivariada.

Uma análise exploratória inicial dos dados foi conduzida utilizando gráficos e medidas descritivas, os quais detectaram padrão espacial das variáveis químicas. Do ponto de vista geoestatístico, o pH e a saturação por bases são dois campos aleatórios gaussianos e latentes. A intuição inicial é que seja ajustado individualmente o modelo (1) para cada vetor de valores observados, que foram notados por $Y_{1}$ e $Y_{2}$ para a saturação por base e o pH, respectivamente. Então $Y_{i}$ segue uma distribuição gaussiana $67$-variada, com matriz de covariâncias $\Sigma_{i}$ e vetor de médias $\mu_{i} = X_{i}\beta_{i}^{*}$, onde $X_{i}$ é uma matriz $67\times p$ com o intercepto e $p-1$ possíveis covariáveis e $\beta_{i}^{*}$ um vetor $p\times1$ com os parâmetros associados a $\mu_{i}$.

Utilizando os resultados da análise exploratória, foram consideradas três formas diferentes para $X_{i}$, média constante, média com tendência induzida pela área de manejo e média com tendência induzida pela coordenada oeste-leste das localizações amostradas. Foi utilizada a família de funções de correlação de Mátern para estruturar a matriz $67\times67$ de correlações, considerando o parâmetro de suavidade $\kappa$ fixo nos valores 0,5, 1, 1.5, 2, 2.5, logo, $\phi_{i}^{*} = \phi_{i}$. Combinando as diferentes escolhas de $\kappa$ e $X_{i}$ foram ajustados diversos modelos. As estimativas paramétricas foram calculadas por máxima verossimilhança, onde foi utilizado a reparametrizaçao $\nu_{i}=\tau_{i}/\sigma_{i}$ e o método numérico de Nelder e Mead (1965) para encontrar as estimativas $\widehat{\phi_{i}}$ e $\widehat{\nu_{i}}$ que maximizam a função de máxima verossimilhança concentrada a partir das quais são encontradas as demais estimativas de parâmetros  de interesse.

Os valores dos máximos estimados da função de máxima verossimilhança concentrada e o Critério de Informação de Akaike (AIC) foram utilizados para selecionar o modelo final de cada atributo. Com os parâmetros estimados dos modelos finais foram conduzidas as krigagens ordinárias para as variáveis químicas do solo. Por último, uma análise de resíduos foi realizada para verificação de pressupostos dos ruídos brancos, os quais se comportaram conforme o esperado.

Após estudo individual de cada variável química, foram ajustados modelos bivariados para o vetor $Y=(Y_{1},Y_{2})$, utilizando as proposições de Diggle e Ribeiro Jr. (2006) e de Matheron (1982). Nas duas abordagens foi utilizado $\mu=X\beta^{*}$, onde $X$ é uma matriz $134\times p$ contento um intercepto para cada resposta e $p-2$ possíveis covariáveis e $\beta^{*}$ é um vetor $p\times1$ com os parâmetros de média associados a $Y$. Para estruturar $\Sigma_{Y}$ somente funções de correlação da família de Matérn foram utilizadas, sendo que o parâmetro $\kappa$ foi sempre considerado fixo em valores similares aos modelos univariados.

Diversos modelos foram ajustados para combinações de diferentes de escolhas de $\kappa$ e $X$. 
Em todos os casos o vetor de parâmetros $\theta$ foi estimado por máxima verossimilhança.
 No BGCCM utilizou-se as reparametrizações $\sigma=\sigma_{01}$,  $\eta=\sigma_{02}/\sigma_{01}$, $\nu_{1}=\sigma_{1}/\sigma_{01}$ e $\nu_{2}=\sigma_{2}/\sigma_{01}$ e no BCRM das reparametrizações $\sigma=\sigma_{1}$, $\nu_{1}=\sigma_{12}/\sigma_{11}$ e $\nu_{2}=\sigma_{22}/\sigma_{11}$. Nas duas abordagens existe forma analítica para os estimadores de $\beta^{*}$ e $\sigma^{2}$, que são função dos demais parâmetros, os quais foram estimados pelo método de Nelder e Mead (1965). Utilizando a propriedade de invariância dos estimadores foi possível encontrar estimativas para todos os parâmetros dos modelos.
Quanto a seleção de modelos e a krigagem, as mesmas ferramentas dos modelos univariados foram utilizadas.

A adoção de modelos geoestatísticos bivariados é uma possível estratégia de análise para esse conjunto de dados, uma vez que, os dois atributos químicos em questão são altamente correlacionados e a coleta da informação sobre a saturação por bases é mais dispendiosa. Como a fazenda continuará sendo monitorada periodicamente, é interessante estruturar a correlação espacial entre e dentre os atributos, pois nas próximas análises de solo é possível amostrar um menor número de pontos georeferenciados para observar a saturação por bases e utilizar a informação do pH para fazer boas inferências. Como exemplo, a amostra da saturação por bases foi dividida em uma subamostra de modelagem e uma subamostra de controle, que contêm 47 e 20 observações, respectivamente. O grupo de modelagem foi utilizado para modelar novamente a saturação por bases sob as abordagens uni e bivariadas. Cabe ressaltar que nos modelos bivariados foi utilizada a informação completa do pH. Com os parâmetros estimados foram calculadas as predições da saturação por bases nas localizações do grupo de controle. Para comparar as metodologias foram calculados os erros de krigagem produzidos em cada caso.

}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "TESALQ"
%%% TeX-command-default: "PdfLaTeX"
%%% End:
