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

\section{Revisão de Literatura}
\subsection{Campos aleatórios}
Um campo aleatório é um processo estocástico que existe em algum espaço real $d$-dimensional, geralmente bi ou
tri-dimensional sua definição é dada por:
\begin{eqnarray*}
\{Z(s_{i}):s_{i}\in G \subset R^{d}\},
\end{eqnarray*}
sendo que $Z(s_{i})$ denota a variável aleatória $Z$ na
localização $s_{i}$ do espaço sob estudo $G$.

Segundo Schmidt e Sansó (2006) e Le e Zidek (2006), a descrição de
um campo aleatório é obtida através das distribuições acumuladas
finito-dimensionais $F$, para qualquer conjunto de localizações $(s_{1}$,$s_{2}$,...,$s_{n})$ pertencentes à região $G$ e
qualquer inteiro $n$:
\begin{eqnarray*}
F_{S_{1},S_{2},...,S_{n}}(z_{1},z_{2},...,z_{n}) \equiv
P(Z(s_{1})\leq z_{1},Z(s_{2})\leq z_{2},...,Z(s_{n})\leq z_{n})
\end{eqnarray*}

Devido a sua simplicidade inferencial, a distribuição de probabilidade gaussiana é uma das mais utilizadas na literatura. Sendo assim, um campo aleatório é dito ser gaussiano se $Z$ segue uma distribuição normal em todas localizações do espaço sob estudo $G$. Para qualquer conjunto finito de localizações $s=(s_{1},s_{2},...,s_{n})$ pertencente a $G$, $Z(s)$ segue uma distribuição normal $n$-variada e é completamente especificado pelo vetor de média $n\times1$, notado por $\mu$, e pela matriz
de covariâncias $n\times n$, notada por $\Sigma$, que no contexto de geoestatística, em geral possui o comportamento de que quanto maior a distância euclidiana entre duas localizações $s_{l}$ e $s_{k}$ quaisquer, menor a correlação entre $Z(s_{l})$ e $Z(s_{k})$. 
A especificação de um modelo com tal comportamento deve
assegurar que a matriz de covariâncias fique positiva definida. 
Diggle e Ribeiro (2007) mostram algumas possíveis construções para obtenção de 
campos aleatórios gaussianos multivariados.
Em geral não é possível ter mais de uma
realização do processo, sendo assim, outras suposições devem ser impostas sobre o campo aleatório gaussiano para a realização de inferências. A restrição mais utilizada é a de que o processo estocástico é estacionário, ou seja, a distribuição de probabilidade associada ao campo aleatório não depende da grandeza de escala das coordenadas, logo, a distribuição conjunta de $(Z(s_{1}), Z(s_{2}),..., Z(s_{n}))$ é igual a distribuição conjunta
de $(Z(s_{1}+h), Z(s_{2}+h),..., Z(s_{n}+h))$, para qualquer incremento $h$. Outra definição menos restritiva é que a média do campo aleatório é constante em toda a região sob estudo e a correlação entre $Z(s_{l})$ e $Z(s_{k})$, para quaisquer $s_{l}$ e $s_{k}$, só depende da distância entre as localizações, ou seja, a grandeza de escala de $Z$ não influencia na estrutura de correlação espacial. Esse tipo de
estacionariedade é conhecido na literatura como estacionariedade
fraca ou de segunda ordem. Uma observação importante é que a
primeira restrição implica na segunda, no entanto, 
o contrário não é válido, exceto no caso do processo espacial ser gaussiano.
No entanto, nem sempre é
fácil verificar as restrições de estacionariedade forte ou fraca,
logo, outra possibilidade menos restritiva é assumir que os
incrementos $[Z(s)-Z(s+h)]$ possuem estacionariedade. Esta
característica é denominada de estacionariedade intrínseca
(SCHANBENBERGER; GOTWAY, 2005). Sendo assim, um campo aleatório é
dito ser intrinsecamente estacionário se para todo $s_{i}$ pertencente a $G$, $E[Z(s_{i})]$=$\mu$ e
$Var[Z(s_{i})-Z(s_{i}+h)]=2\gamma(h)$, sendo $\gamma(h) = Var(Z(s_{i})) - Cov(Z(s_{i});Z(s_{i}+h))$ e 
denominado de semivariograma.

Em alguns casos a suposição de estacionariedade não é válida, sendo assim, diversas abordagens são propostas para contornar esse problema. Quando a média do processo estocástico não é constante na região sob estudo, geralmente, utiliza-se a inclusão de covariáveis na modelagem, onde a média é tratada como efeito fixo, e sua interpretabilidade é igual a de modelos lineares. 
Já para problemas com variâncias e covariâncias não constantes, uma técnica mais simples é a utilização de transformação nos valores observados do campo aleatório (CHRISTENSEN; DIGGLE; RIBEIRO Jr., 2001).
Alternativas incluem modelos não estacionários que, embora mais gerais, requerem um grande volume de informação
para serem identificados e 
dificilmente podem ser
estimados para conjuntos esparsos de dados.

Mesmo que um campo aleatório possua algum tipo de estacionariedade, o padrão espacial de correlação pode depender das distâncias e das direções envolvidas entre as localizações e 
a essa característica é dado o nome de anisotropia. 
Em problemas práticos de geoestatística não é fácil identificar tal característica 
a partir dos dados observacionais, No entanto, informações contextuais sobre 
a natureza de algumas variáveis podem indicar a utilização de tal abordagem, como por exemplo, estudos sobre poluição, onde o padrão dos ventos gera correlação espacial dependente das direções entre localizações amostradas. Silva (2006) mostra algumas técnicas utilizadas em problemas de anisotropia geométrica.

Um campo aleatório gaussiano é dito ser homogêneo se ele for estacionário e seu padrão de correlações não depende das direções, ou seja, isotrópico. 
Utilizando essa suposição, o processo estocástico, embora restritivo, 
consegue modelar diversos problemas naturais.

Com a suposição de homogeneidade de um campo aleatório gaussiano, é necessário estabelecer uma função matemática que dependa apenas das distâncias entre as localizações amostradas do espaço sob estudo e que estruture a matriz de covariâncias, de forma que esta seja positiva definida e com o comportamento típicamente 
utilizado em geoestatística. 
Encontrar tais funções que sejam válidas pode ser um procedimento 
complexo. Na prática, o usual é adotar famílias de funções ditas permissíveis, 
entre propostas na literatura.
A seguir são descritas duas famílias que têm como característica possuirem dois parâmetros, um deles controlando
a extensão da dependência espacial e o outro controlando a suavidade do processo.

\textbf{Função de Matérn}

Essa família de funções de correlação foi proposta por Bertil Matèrn (1986) e
possui a seguinte forma:
\begin{eqnarray*}
\rho(h)=\left{2^{\kappa-1}\Gamma(\kappa)\right}^{-1}(h/\phi)^{\kappa}K_{\kappa}(h/\phi),
\end{eqnarray*}
sendo $\Gamma(\cdot)$ a função gama, $K_{\kappa}(\cdot)$ é a função Bessel de ordem $\kappa$,  
$h$ a distância euclidiana entre duas localizações quaisquer do campo aleatório, 
os parâmetros dessa função são $\phi$ > 0 e $\kappa$ > 0.
O primeiro desses parâmetros está relacionado ao alcance das correlações, sendo que maiores
valores indicam dependências espaciais de maior alcance.
O segundo parâmetro é relacionado à suavidade
do processo, de forma que quanto maior $\kappa$ maior a suavidade.

\textbf{Função Exponencial Potência}
Essa família de funções de correlação 
possui características similares às da função de Matèrn
onde $\phi$ e $\kappa$ são interpretados de forma análoga.
Entretanto, neste caso, $\kappa$ é limitado ao intervalo $[0,2]$ 
e não definindo uma sequência de diferenciabilidade dos processos associados que passam de
não diferenciáveis para $\kappa < 2$ para infinitamente diferenciável  para $\kappa > 2$.
A forma desta funçõa de correlação é:
\begin{eqnarray*}
\rho(h)=\exp{(h/\phi)^{\kappa}},
\end{eqnarray*}

Essas duas funções são muito utilizadas devido a capacidade de produzir comportamentos distintos quanto a suavidade do processo, ou seja, é possível modelar processos mais ou menos diferenciáveis. Além disso, sob estacionariedade fraca, essas funções possuem propriedade conhecidas e desejáveis.
Schabenberger e Gotway (2005) apresentam e discutem detalhes sobre tais propriedades. 
Silva (2006), em um trabalho discutindo a construção de funções de correlação válidas para 
modelos fgeoestatísticos espaço temporais 
 apresenta diversas outras funções de correlação conhecidamente válidas e bem como uma extensa lista de
referências sobre o tema.

\subsection{Modelos geoestatísticos gaussianos univariados}
Considerando que em alguma área $G$ exista um campo aleatório gaussiano $Z$ latente, ou seja, o processo existe mas não é observável, é necessário fazer uma amostragem de $n$ localizações espaciais dentro da área $G$ e observar valores de atributo de interesse nas localizações amostradas. Sendo assim, existe um vetor $Y(s)$  de 
$n\times1$ valores observados em $s=(s_{1},s_{2},...,s_{n})$, que segundo Diggle e Ribeiro Jr (2007),
sob o prassuposto de gaussianidade, pode ser modelado como:
\begin{eqnarray}
 Y(s) = \mu + Z(s) + \epsilon,
\end{eqnarray}
onde $\mu=X\beta$, onde $X$ é uma matriz $n\times q$ contendo $q-1$ possíveis covariáveis, $\beta$ um vetor $q\times1$ de parâmetros associados a $X$, $Z(s)$ um campo aleatório gaussiano que possui vetor de médias $n\times1$ nulo e matriz de covariâncias $\Sigma$ de dimensão $n\times n$, onde cada elemento $\Sigma_{i,j}$ é igual a $Cov(Z(s_{i});Z(s_{j}))$, para todo $s_{i}$ e $s_{j}$ pertencentes a $s$, e $\epsilon$ um vetor $n$x$1$ de ruídos brancos, que por suposição são independentes e identicamente distribuidos com distribuição de probabilidade normal com média zero e desvio padrão $\tau$.

Utilizando (1) é possível encontrar a distribuição de probabilidade de $Y(s)$, que é gaussiana $n$-variada, com vetor de médias $X\beta$ e matriz de covariâncias $\Sigma_{Y} = \Sigma + \tau^{2} I$, onde $I$ é uma matriz identidade $n\times n$. Logo, existe um vetor de parâmetros $\theta=(\beta,\sigma^{2},\phi^{*},\tau^{2})$ a ser estimado, onde $\phi^{*}$ é um vetor de parâmetros associados a função de correlação utilizada. Cabe observar que $Cov(Z(s_{i});Z(s_{j}))=\sigma^{2} \rho(Z(s_{i});Z(s_{j}))$, para todo $s_{i}$ e $s_{j}$ pertencentes a $s$ e $\sigma^{2}$>0.
\subsubsection{Estimação dos parâmetros}
Especificado um modelo adequado para os dados em questão, o próximo passo é fazer a
estimação dos parâmetros. Se o campo aleatório é intrinsecamente
estacionário pode-se trabalhar com uma estimação para o semivariograma.
A  expressão de uma estimativa empírica para o
semivariograma através dos estimadores do método de momentos é:
\begin{eqnarray}
\widehat{\gamma}(h)= \frac{\sum_{|N(h)|}(Z(s_{i})-Z(s_{j}))^{2}}{2|N(h)|}
\end{eqnarray}
em que $|N(h)|$ é o número de pontos abrangidos pela
distância $h$.

Devido a relação entre o semivariograma empírico e as funções de
correlação válidas, muitos trabalhos aplicados de geoestatística
utilizam um modelo em função dos parâmetros de variabilidade e correlação
que se ajuste aos valores calculados em (2), isso pode
ser feito por meio de métodos que ajustam um modelo de função de semivariância aos dados, sendo que 
os parâmetros de tal função correspondem aos parâmetros do modelo. 
Essa abordagem para
estimar os parâmetros de variabilidade e da função de correlação podem não ser muito
precisas, pois os valores dos semivariogramas empíricos podem se
afastar muito do semivariograma real e desconhecido, devido ao
tamanho, acaso amostral, variabilidade e subjetividade do variograma empírico e a intrínsica correlação entre seus pontos. Diggle e Ribeiro Jr. (2007) discutem com mais detalhes esse assunto e recomendam que 
esse tipo de estimação deve ser utilizado, 
na maioria dos casos, apenas como análise descritiva inicial, 

Por outro lado, assumindo que o campo aleatório possui
estacionariedade forte, pode-se optar por estimadores baseados na
verossimilhança ou verossimilhança restrita. Esta abordagem  
consiste em utilizar os valores observados da variável resposta para encontrar um vetor $\widehat{\theta}$ que seja o ponto de máximo da função de verossimilhança associada a $\theta$. 
Por simplicidade matemática, normalmente  utiliza-se o logaritmo da função de verossimilhança para fazer a estimação, que associada a (1) tem a seguinte forma:
\begin{eqnarray*}
l(\theta;Y(s))=-0.5(n\ln(2\pi)+\ln(|\Sigma_{Y}|)+(Y(s)-X\beta)^{t} \Sigma_{Y}^{-1} (Y(s)-X\beta)).
\end{eqnarray*}

Maiores detalhes sobre técnicas e propriedades da estimação por máxima verossimilhança são expostos por Azzalini (1996) e Bickel e Doksum (1976). No contexto de geoestatística, Diggle e Ribeiro Jr. (2007) adotam a utilização da reparametrização $\nu=\tau / \sigma$, a qual facilita a estimação de $\theta$. O vetor de parâmetros ser estimado passar a ser $\theta^{*}=(\beta,\sigma^{2},\phi^{*},\nu^{2})$ e $\Sigma_{Y}$ pode ser escrita como $\sigma^{2}V$, onde $V$ é uma matriz $n\times n$ que depende apenas de $\nu$ e $\phi^{*}$. O logaritmo da função de verossimilhança fica da seguinte forma:
\begin{eqnarray}
l(\theta^{*};Y(s)) \propto -0.5(n \ln(\sigma^{2}) - \ln(|V|) - \sigma^{-2} (Y-X\beta)^{t} V^{-1} (Y-X\beta)
\end{eqnarray}
sendo que existem formas analíticas apenas para os estimadores de $\beta$ e $\sigma^{2}$:
\begin{eqnarray*}
\widehat{\beta} = (X^{t}V^{-1}X)^{-1} (X^{t}V^{-1}Y) \; ; \; \widehat{\sigma^{2}} = n^{-1} (Y-X\beta)^{t} V^{-1} (Y-X\beta)
\end{eqnarray*}

As expressões de $\widehat{\beta}$ e $\widehat{\sigma^{2}}$ são funções dos demais parâmetros e, 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 $(\widehat{\beta},\widehat{\sigma^{2}})$ (DIGGLE; RIBEIRO Jr., 2007).
Para $\phi^{*}$ e $\nu^{2}$ não existem formas analíticas para os estimadores. Utilizando $\widehat{\beta}$ e $\widehat{\sigma^{2}}$ em (3), tem-se o logaritmo da função de verossimilhança concentrada, que depende apenas de $\theta_{c}=(\nu,\phi^{*})$.

Para encontrar $\widehat{\theta_{c}}$ é possível utilizar métodos numéricos de maximização de funções, como por exemplo, o método de Nelder e Mead (1965), o qual calcula numericamente $\widehat{\phi^{*}}$ e $\widehat{\nu^{2}}$ e a matriz Hessiana estimada, denotada por $H$. Com os parâmetros da função de máxima verossimilhança concentrada estimados, por invariância é possível encontrar as estimativas de $\beta$ e $\sigma^{2}$ e $\tau^{2}$.

Utilizando as propriedades assintóticas dos estimadores de máxima verossimilhança e o método delta pode-se encontrar que a distribuição de probabilidade de $\widehat{\theta}$ é $N( \theta ; \Sigma_{\theta} )$, sendo $\Sigma_{\theta} = \Delta^{t} \Sigma_{\theta^{*}} \Delta$, onde a $i$-ésima coluna de $\Delta$ é o vetor $\frac{\partial l(\theta_{i})}{\partial \theta^{*}}$ e:
\[ \Sigma_{\theta^{*}} = \left[\begin{array}{cc}
\Sigma_{\beta,\sigma^{2}}   & O     \\
O^{t} & \Sigma_{\theta_{c}} \end{array} \right] \]
\noindent sendo $\Sigma_{\beta,\sigma^{2}}$ a matriz de covariâncias de $(\widehat{\beta},\widehat{\sigma^{2}})$, que possui forma analítica, $\Sigma_{\theta_{c}}= -H^{-1}$ é a matriz de covariâncias de $\widehat{\theta_{c}}$, e $O$ uma matriz de zeros, pois não é possível recuperar as covariâncias entre $(\hat{\beta},\hat{\sigma}^{2})$ e $\hat{\theta}_{c}$.

Geralmente o objetivo final dos estudos com dados georeferenciados é calcular as predições espaciais em localizações não amostradas. Esse processo de predição é denominado de krigagem e é baseado nas propriedades da distribuição normal multivariada.

\subsubsection{Krigagem}
A krigagem nada mais é do que o processo de predição do campo aleatório em localizações não amostradas. O nome krigagem é uma homenagem ao pesquisador sul-africano D.G. Krige que foi um dos pioneiros em estudos de predição espacial em inas de ouro da África do Sul.

Suponha que o interesse é fazer predição para $Z$ nas localizações $s^{*}$, ou seja, devemos fazer predição para $Z(s^{*})$ tal que o erro quadrático médio seja mínimo. A distribuição de probabilidade de $(Z(s^{*})|Y(s))$ gera as predições com melhor precisão. Em particular a esperança desta distribuição minimiza o erro quadrático médio de predição e a variância fornece a precisão da estimativas. Sendo assim pode-se obter:
\begin{eqnarray}
E(Z(s^{*})|Y(s)) = \mu_{Z} + \Sigma_{Z,Y} \Sigma_{Y}^{-1} (Y(s)-X\beta)
\end{eqnarray}
sendo $\mu_{Z}$ a média de $Z(s^{*})$ e $\Sigma_{Z,Y}$ é a matriz de covariâncias cruzadas entre $Z(s^{*})$ e $Y(s)$.
Além disso a variância preditiva é conhecida:
\begin{eqnarray}
Var(Z(s^{*})|Y(s)) = \Sigma_{Z} - \Sigma_{Z,Y} \Sigma_{Y}^{-1} \Sigma_{Z,Y}^{t}
\end{eqnarray}
sendo $\Sigma_{Z}$ a matriz de covariâncias de $Z(s^{*})$.
Cabe ressaltar que as expressões (4) e (5) são obtidas através das propriedades da distribuição gaussiana multivariada.

Do ponto de vista geoestatístico, as krigagens mais frequentemente utilizadas são a simples 
e a ordinária, que utilizam (3) e (4), mas que se diferenciam quanto a suposição de conhecimento prévio
sobre os parâmetros da média do processo (DIGGLE; RIBEIRO Jr., 2007).



\section{Modelos geoestatísticos gaussianos bivariados}
Quando existem dois campos aleatórios gaussianos de interesse, a intuição inicial é que seja ajustado um modelo geoestatístico univariado para cada vetor de valores observados. Porém, em alguns casos existe correlação estatística entre as variáveis aleatórias, o que leva a possibilidade de adotar modelos geoestatísticos bivariados. Cabe ressaltar que somente com essa justificativa estatística não é vantajoso utilizar tal abordagem, é necessário que existam justificativas e vantagens práticas para esse aumento de complexidade dos modelos geoestatísticos, que nesses casos devem estruturar as correlações entre e dentre os campos aleatórios.

Supondo a existência de dois campos aleatórios gaussianos de interesse, pode-se utilizar seguinte o modelo:
\begin{eqnarray}
Y_{i} = \mu_{i} + Z_{i}, \ i=1,2
\end{eqnarray}
sendo $Y_{i}$ um vetor $n_{i}\times1$ de valores observados do campo aleatótio gaussiano latente $Z_{i}$, que possui vetor de médias nulo $n_{i}\times1$ e matriz de covariâncias $\Sigma_{i}$ $n_{i}\times n_{i}$, $\mu_{i}$ é um vetor, que possui os parâmetros de médias associados a $Y_{i}$. Observe que em (6) não existem termos de ruídos brancos.

A distribuição de $Y=(Y_{1},Y_{2})$ é de interesse final e possui distribuição gaussiana $n$-variada, sendo $n=n_{1}+n_{2}$, com vetor de médias $\mu=(\mu_{1},\mu_{2})$ e matriz de covariâncias $\Sigma_{Y}$, positiva definida que possui o comportamento empírico de correlações utilizado em geoestatística e que pode ser particionada da forma:
\[ \Sigma_{Y} = \left[\begin{array}{cc}
\Sigma_{1}   & \Sigma_{1,2}     \\
\Sigma_{1,2}^{t} & \Sigma_{2} \end{array} \right] \]
\noindent onde $\Sigma_{i}$ é uma matriz $n_{i}\times n_{i}$ das covariâncias da variável $Y_{i}$, $i=1,2$, e $\Sigma_{1,2}$ uma matriz $n_{1}\times n_{2}$ com as covariâncias cruzadas entre $Y_{1}$ e $Y_{2}$.

Nesse contexto a especificação do modelo exige  
uma matriz $\Sigma_{Y}$ válida. 
A abordagem por contrução via modelos separáveis é a mais simples para resolver 
esse problema e utiliza decomposições os termos latentes de (6).

\subsection{Modelo gaussiano bivariado com componente de correlação parcialmente comum}
Diggle e Ribeiro Jr. (2006) propuseram uma abordagem para problemas geoestatísticos bivariados, a qual é abreviada nesse trabalho por BGCCM (\textit{bivariate Gaussian common component model}
e que utiliza as seguintes decomposições dos campos aleatórios latentes de (6):
\begin{eqnarray}
Z_{i}=\sigma_{0i}S_{0}+\sigma{i}S_{i}, \ i=1,2,
\end{eqnarray}
sendo $\sigma^{*}=(\sigma_{01},\sigma_{1},\sigma_{02},\sigma_{2})$ um vetor de parâmetros de dispersão associados a (6) e $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. 
As correlações destes termos, quando combinadas desta forma, vão gerar correlações cruzadas válidas 
entre as variáveis observáveis 
$Y_{1}$ e $Y_{2}$, uma vez que, o campo aleatório $S_{0}$ é comum as duas respostas.
Portanto, sob este modelo,
três funções de correlação válidas para $S_0$, $S_1$ e $S_2$ devem ser adotadas 
para gerar uma matriz $\Sigma_{Y}$ válida.

Logo, o modelo dado por (6) fica da seguinte maneira:
\[ \left\{ \begin{array}{l}
Y_{1}=\mu_{1}+\sigma_{01}S_{0}+\sigma_{1}S_{1} \\
Y_{2}=\mu_{2}+\sigma_{02}S_{0}+\sigma_{2}S_{2} \end{array}\right.\]

Definindo $Y_{i}(s_{l})$ e $Y_{i}(s_{k})$ como observações de $Y_{i}$ feitas em duas localizações quaisquer $s_{l}$ e $s_{k}$, que estão separados por uma distância euclidiana $h$, para todo $l,k=1,2,\ldots,n_{i}$ e $i=1,2$, tem-se que o elemento $\Sigma_{i,(l,k)}$ é dado por $Cov(h)=\sigma_{0i}^{2}\rho_{0}(h)+\sigma_{i}^{2}\rho_{i}(h)$, sendo $\rho_{0}$ e $\rho_{i}$ as funções de correlação adotadas para $S_{i}, i=0,1,2$, respectivamente. Utilizando propriedades básicas de covariâncias, pode-se encontrar $\Sigma_{1,2}$, que é igual a $\sigma_{01}\sigma_{02}R_{0}$, sendo $R_{0}$ uma matriz $n_{1}\times n_{2}$ das correlações cruzadas entre as respostas e depende da função de correlação adotada para $S_{0}$. Dessa forma $\Sigma_{Y}$ fica completamente estruturada. Por simplicidade, suponha que existam apenas duas localizações amostradas $s_k$ e $s_l$ sendo $h = h_{kl}$ a distância entre elas. As covariâncias para um vetor $(Y_i(s_l), Y_i(s_k), Y_j(s_l), Y_j(s_k))\prime$ são dadas por:
\begin{center}
\[ \Sigma_{Y} = \left[\begin{array}{cccc}
\sigma_{01}^{2}+\sigma_{1}^{2} & \sigma_{01}^{2}\rho_{0}(h) +\sigma_{1}^{2}\rho_{1}(h) & \sigma_{01}\sigma_{02} & \sigma_{01}\sigma_{02}\rho_{0}(h) \\
\sigma_{01}^{2}\rho_{0}(h) +\sigma_{1}^{2}\rho_{1}(h) & \sigma_{01}^{2}+\sigma_{1}^{2} & \sigma_{01}\sigma_{02}\rho_{0}(h) & \sigma_{01}\sigma_{02} \\
\sigma_{01}\sigma_{02} & \sigma_{01}\sigma_{02}\rho_{0}(h) & \sigma_{02}^{2}+\sigma_{2}^{2} & \sigma_{02}^{2}\rho_{0}(h) +\sigma_{2}^{2}\rho_{2}(h) \\
\sigma_{01}\sigma_{02}\rho_{0}(h) & \sigma_{01}\sigma_{02} & \sigma_{02}^{2}\rho_{0}(h) +\sigma_{2}^{2}\rho_{2}(h)  & \sigma_{02}^{2}+\sigma_{2}^{2}  \end{array} \right] \].
\end{center}

\ \

Sendo assim, a distribuição de probabilidade do vetor $Y=(Y_{1},Y_{2})$ está estabelecida e depende do vetor de parâmetros $\theta=(\beta^{*},\sigma^{*},\phi_{0}^{*},\phi_{1}^{*},\phi_{2}^{*})$, onde $\beta^{*}$ é um vetor de parâmetros associado a $\mu$ e $\phi_{j}^{*}$ é um vetor de parâmetros associado a escolha da função $\rho_{j}, j=0,1,2$. Para fazer a estimação dos parâmetros e a predição espacial pode-se utilizar as mesmas técnicas envolvidas nos modelos geoestatísticos univariados. No entanto, é necessário utilizar reparametrizações adicionais em $\theta$ para encontrar a função de verossimilhança concentrada.

\subsection{Modelo bivariado de co-regionalização}
Gelfand et al. (2004) propôs uma abordagem para encontrar uma estrutura paramétrica válida para $\Sigma_{Y}$, que é abreviada aqui por BCRM (\textit{bivariate coregionalisation model}). 
Nesta abordagem os termos latentes de (6) são decompostos da seguinte forma:
\[ \left\{ \begin{array}{l}
Z_{1}=\sigma_{11}S_{1} \\
Z_{2}=\sigma_{12}S_{1} + \sigma_{22}S_{2} \end{array}\right.\]
\noindent sendo $S_{1}$ e $S_{2}$ campos aleatórios gaussianos mutuamente independentes, com vetores de médias nulo, variância unitária e com correlações determinadas pela escolha de funções de correlação conhecidamente válidas. Como campo aleatório $S_{1}$ é comum as duas variáveis, é induzida a correlação empirica entre as respostas. Gelfand et al. (2004) mostram mais detalhes sobre essa abordagem do ponto de vista de modelos geoestatísticos.

Utilizando o BCRM a expressão (6) do modelo para as observasções  é representada da seguinte forma:
\[ \left\{ \begin{array}{l}
Y_{1}=\mu_{1}+\sigma_{11}S_{1}\\
Y_{2}=\mu_{2}+\sigma_{12}S_{1}+\sigma_{22}S_{2} \end{array}\right.\]

Definindo $Y_{i}(s_{l})$ e $Y_{i}(s_{k})$ como observações feitas da resposta $Y_{i}$ em duas localizações quaisquer $s_{l}$ e $s_{k}$, separadas por uma distância euclidiana $h = h_{lk}, l,k=1,2,...,n_{i}$ e $i=1,2$, tem-se que o elemento $\Sigma_{1,(l,k)}$ é dado por $Cov(h)=\sigma_{11}^{2}\rho_{1}(h)$ e o elemento $\Sigma_{2,(l,k)}$ é dado por $Cov(h)=\sigma_{12}^{2}\rho_{1}(h) + \sigma_{22}^{2}\rho_{2}(h)$, onde $\rho_{1}$ e $\rho_{2}$ são as funções de correlação adotadas para $S_{1}$ e $S_{2}$, respectivamente. Utilizando propriedades básicas de covariâncias, pode-se encontrar $\Sigma_{1,2}$, que é igual a $\sigma_{11}\sigma_{12}R_{1}$, onde $R_{1}$ é uma matriz $n_{1}\times n_{2}$ das correlações cruzadas entre as respostas e depende apenas da função de correlação adotada para $S_{1}$. Dessa forma $\Sigma_{Y}$ fica completamente estruturada. Por simplicidade, suponha que existam apenas duas localizações amostradas, logo o vetor  $(Y_i(s_l), Y_i(s_k), Y_j(s_l), Y_j(s_k))\prime$ tem matrix de covariância:

\begin{center}
\[ \Sigma_{Z} = \left[\begin{array}{cccc}
\sigma_{11}^{2} & \sigma_{11}^{2}\rho_{1}(h) & \sigma_{11}\sigma_{12} & \sigma_{11}\sigma_{12}\rho_{1}(h) \\
\sigma_{11}^{2}\rho_{1}(h) & \sigma_{11}^{2} & \sigma_{11}\sigma_{12}\rho_{1}(h) & \sigma_{11}\sigma_{12} \\
\sigma_{11}\sigma_{12} & \sigma_{11}\sigma_{12}\rho_{1}(h) & \sigma_{12}^{2}+\sigma_{22}^{2} & \sigma_{12}^{2}\rho_{1}(h) +\sigma_{22}^{2}\rho_{2}(h) \\
\sigma_{11}\sigma_{12}\rho_{1}(h) & \sigma_{11}\sigma_{12} & \sigma_{12}^{2}\rho_{1}(h) +\sigma_{22}^{2}\rho_{2}(h)  & \sigma_{12}^{2}+\sigma_{22}^{2}  \end{array} \right] \].
\end{center}
\ \

Sendo assim, a distribuição de probabilidade do vetor $Y=(Y_{1},Y_{2})$ está estabelecida e depende do vetor de parâmetros $\theta=(\beta^{*},\sigma^{*},\phi_{1}^{*},\phi_{2}^{*})$, onde $\beta^{*}$ é um vetor de parâmetros associado a $\mu$ e $\phi_{j}^{*}$ é um vetor de parâmetros associado a escolha da função $\rho_{j}$, para todo $j=1,2$. Para fazer a estimação dos parâmetros e a predição espacial pode-se utilizar as mesmas técnicas envolvidas nos modelos geoestatísticos univariados com reparametrizações adequadas em $\theta$ para encontrar a função de verossimilhança concentrada.

Uma generalização deste modelo consiste na adição de um termo de ruídos não correlacionados
que podem ser associados a todos ou alguns dos $Y_i$. Neste caso a equação básica dos modelos dada em (6)
fica extendida na forma a seguir onde $U$ é um processo ruído branco,
\begin{eqnarray}
Y_{i} = \mu_{i} + Z_{i} + U_{i}, \ i=1,2.
\end{eqnarray}


\subsection{Inferência e predição Bayesiana}
%No contexto de modelagem estatística, a estimação os parâmetros desconhecidos é uma fase importante das pesquisas, para executar tal tecnicalidade pode-se utilizar diversas ferramentas estatísticas, como por exemplo, estimadores de máxima verossimilhança. No entanto, em problemas com estruturas paramétricas mais complexas, as técnicas frequentistas não são facilmente aplicáveis, sendo assim, uma abordagem concorrente é a utilização de técnicas bayesianas.

No contexto de geoestatística os parâmetros dos modelos são estimados e utilizados na obtenção de predições 
espaciais. Existe portanto incerteza sobre estes parâmetros que podem ser incorporadas na 
incerteza associada às predições espaciais.
O enfoque bayesiano fornece um ferramental natural para tal, representando a incerteza sobre os parâmetros
através de distribuições de probabilidades.

A metodologia de estimação por máxima verossimilhança utiliza apenas as informações provinientes dos dados, do modelo de probabilidade adotado e das teorias assintóticas para fazer inferências. 
Já no contexto bayesiano também é utilizado a informação que existe a $\emph{priori}$ sobre os parâmetros, ou seja, os parâmetros são desconhecidos, mas pode existir um conhecimento inicial sobre os mesmos
expressos por meior das distribuições a priori, ou seja, na prática 
não existe diferenciação entre parâmetros e observações, sendo ambos tratados como quantidades aleatórias.

A definição da informação a $\emph{priori}$ dos parâmetros é fonte de debate entre os estatísticos, uma vez que, é considerada uma determinada distribuição de probabilidade aos parâmetros e cada pesquisador pode considerar informações iniciais diferentes, dependendo do grau de conhecimento sobre o problema, ou seja, as distribuições de probabilidade associadas aos parâmetros podem ser mais ou menos informativas.
Naturalmente,  diferentes suposições levam a modelos de probabilidades diferentes e consequentemente a resultados inferenciais e predições diferentes.

O teorema de Bayes quantifica a atualização da informação sobre os parâmetros, que é revisada após a observação de uma quantidade aleatória $Y=y$:
\begin{eqnarray*}
p(\theta | y) = \frac{p(\theta , y)}{p(x)} = \frac{p(y | \theta) p(\theta)}{p(y)}
\end{eqnarray*}
sendo $\theta$ o parâmetro ou conjunto de parâmetros de interesse, uma observação importante é que $p(y)$ é uma integral de $p(\theta , y)$ no espaço paramétrico de $\theta$, ou seja, esse termo é uma constante.

Observado o valor $y$, $p(y | \theta)$ fornece uma plausibilidade para cada valor possível de $\theta$, e por isso é chamada de função de verossimilhança do parâmetro de interesse; $p(\theta)$ é a distribuição a $\emph{priori}$ dos parâmetros, que é arbitrária e definida pelo pesquisador. Desta forma temos a definição da distribuição a $\emph{posteriori}$, que nada mais é do que a plausibilidade dos possíveis valores de $\theta$ para um determinado conjunto de dados atualizada, ponderada, pela informação inicial sobre os parâmetros, logo:
\begin{center}
$p(\theta | y) \propto L(\theta ; y) p(\theta)$
\end{center}
\noindent $L(\theta ; y)$ é a função de verossimilhança associada ao modelo sob estudo, e o símbolo $\propto$ indica uma proporção, ou seja, os termos constantes do teorema de Bayes são considerados apenas como constantes normalizadora, tal utilização facilita e não interfere na estimação dos parâmetros. No entanto, 
para outros casos, como seleção de modelo, se faz necessário recuperar os valores das constantes, o que de certa forma é bastante intuitivo, como $p(\theta | y) = \kappa p(y | \theta) p(\theta)$ e se essa equação for integrada com relação a $\theta$, $p(\theta | x) = 1$ e a constante normalizadora $\kappa$ pode ser encontrada:
\begin{center}
$\kappa^{-1} = \int p(y | \theta) p(\theta) d\theta = E_{\theta}[p(Y | \theta] = p(y)$
\end{center}
\noindent chamada de distribuição preditiva, ou seja, é a distribuição esperada para a observação $y$ dado $\theta$.

Em algumas aplicações como as de geoestatística, após observado $Y=y$ existe interesse em fazer previsões para uma quantidade $Y_0$ não observada e  relacionada a $\theta$ da seguinte forma:

\begin{center}
$p(y_0 | y) = \int p(y_0 , \theta | y) d\theta = \int p(y_0 | \theta , y) p(\theta | y) d\theta.$
\end{center}
\noindent na maioria dos casos essa integral não possui solução analítica, logo métodos de aproximação podem ser utilizados.

Após definida a distribuição a $\emph{posteriori}$ dos parâmetros é possível resumir informações sobre os parâmetros através de técnicas inferenciais, utilizando normalmente probabilidades e esperanças, que dependendo da complexidade do problema não são resolvidas analiticamente. Nesta etapa que se faz necessário utilização de métodos computacionais intensivos como, por exemplo, o Monte Carlo via Cadeias de Markov, denominada de MCMC. No contexto de geoestatística o amostrador de Gibbs e o algoritmo Metropolis-Hastings são os métodos mais utilizados, Gamerman (1997) apresenta uma visão geral sobre esses métodos.

Em modelos geoestatísticos Kitanidis (1986), Handcock e Stein (1993), Diggle e Ribeiro (2002), Banerjee et. al. (2004), entre outros adotam  a abordagem Bayesiana em modelos univariados. 
Schmidt e Gelfand (2003), Finley et al. (2007) e Gelfand et. al. (2007) adotam esta abordagem em modelos multivariados.
No caso de modelos 
A incorporação de procedimentos bayesianos é uma extensão natural 
e vem sendo adaptada aos modelos bivariados considerados aqui
como extensão deste trabalho.

}


@article{Finley:Banerjee:Carlin:2007:JSSOBK:v19i04,
  author =	"Andrew O. Finley and Sudipto  Banerjee  and Bradley  P.  Carlin",
  title =	"spBayes: An R Package for Univariate and Multivariate Hierarchical Point-referenced Spatial Models",
  journal =	"Journal of Statistical Software",
  volume =	"19",
  number =	"4",
  pages =	"1--24",
  day =  	"24",
  month =	"4",
  year = 	"2007",
  CODEN =	"JSSOBK",
  ISSN = 	"1548-7660",
  bibdate =	"2007-04-24",
  URL =  	"http://www.jstatsoft.org/v19/i04",
  accepted =	"2007-04-24",
  acknowledgement = "",
  keywords =	"",
  submitted =	"2006-10-29",
}


@Book{banerjee+carlin+gelfand:2004,
  author =	 {Banerjee, S. and Carlin, B. P. and Gelfand, A. E.},
  title =	 {Hierachical Modelling and Analysis for spatial Data},
  publisher =	 {Chapman and Hall},
  address =	 {Boca Raton},
  year =	 {2004},
}


@Article{gelfand+schmidt+banerjee+sirmans:2004,
  author = 	 {Alan E. Gelfand and Alexandra M. Schmidt and Sudipto Banerjee and C. F. Sirmans},
  title = 	 {Nonstationary multivariate process modeling through spatially varying coregionalization (with Discussion)},
  journal = 	 {Test},
  year = 	 {2004},
  OPTkey = 	 {},
  volume = 	 {13},
  OPTnumber = 	 {2},
  pages = 	 {263--312},
  OPTmonth = 	 {},
  OPTnote = 	 {},
  OPTannote = 	 {}
}


@Article{schmidt+gelfand:2003,
  author = 	 {Alexandra M. Schmidt and Alan E. Gelfand},
  title = 	 {A Bayesian corregionalization approach for multivariate pollutant data},
  journal = 	 {Journal of Geophysical Research --- Atmospheres},
  year = 	 {2003},
  OPTkey = 	 {},
  volume = 	 {108 (D24)},
  pages = 	 {8783},
  OPTmonth = 	 {},
  OPTnote = 	 {},
  OPTannote = 	 {}
}

@Article{kitanidis:1986,
  author = 	 {Kitanidis, P. K.},
  title = 	 {Parameter uncertainty in estimation of spatial functions: {B}ayesian analysis},
  journal = 	 WRR,
  year = 	 {1978},
  volume = 	 {22},
  OPTnumber = 	 {4},
  pages = 	 {499--507},
}
@Article{handcock+stein:1993,
  author = 	 {Handcock, M. and Stein, M.},
  title = 	 {A {B}ayesian analysis of kriging},
  journal = 	 {Technometrics},
  year = 	 {1993},
  volume = 	 {35},
  OPTnumber = 	 {},
  pages = 	 {403--410},
}

@Article{diggle+ribeiro:2002,
  author = 	 {Diggle, P. J. and {Ribeiro Jr}, P. J.},
  title = 	 {{B}ayesian inference in {G}aussian model based geostatistics},
  journal = 	 {Geographical and Environmental Modelling},
  year = 	 {2002},
  volume = 	 {6},
  OPTnumber = 	 {2},
  pages = 	 {129--146},
}

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