1ª Escola de Verão - GET/IME/UFF

Parte 3

1 Superfícies Contínuas ou Geoestatística

  • Considera-se que o processo aleatório a ser estudado varia ao longo de um espaço contínuo \(\mathbb{R}^p\), onde \(p = 1,2,3\), usualmente. As localizações observadas do processo são conhecidas e podem ser identificadas pela latitude, pela longitude e pela altitude, por exemplo.

  • Seja \(\mathbb{G}\) um subconjunto de \(\mathbb{R}^p\) com volume positivo \(p\)-dimensional. Seja \(\{H(\boldsymbol{s})\), \(\boldsymbol{s} \in \mathbb{G}\}\), um processo aleatório. Por exemplo, suponha que o fenômeno de interesse \(H(\boldsymbol{s})\) seja a temperatura marcada em um termômetro em uma localização fixa sendo \(\boldsymbol{s} = (s_1, s_2)\) onde \(s_1\) é a latitude e \(s_2\) é a longitude onde estas temperaturas são observadas. Neste caso, tem-se que \(p=2\).

  • Suponha que a média do processo aleatório \(\{H(\boldsymbol{s})\), \(\boldsymbol{s} \in \mathbb{G}\}\) existe e denote-a por \(\mu(\boldsymbol{s}) = E[H(\boldsymbol{s})]\). Suponha que a variância, \(Var[H(\boldsymbol{s})]\), também existe para todo \(\boldsymbol{s} \in \mathbb{G}\).

  • O processo é dito ser Gaussiano quando a distribuição conjunta de \(\boldsymbol{H} = (H(\boldsymbol{s_1}), \ldots, H(\boldsymbol{s_n}))^{'}\) for a distribuição normal multivariada, para qualquer \(n=1, 2, \ldots\) e para quaisquer coleções \(\{\boldsymbol{s_1}, \ldots, \boldsymbol{s_n}\}\), onde \(\boldsymbol{s_i} \in \mathbb{G}\).

  • O processo é dito ser estritamente estacionário se a distribuição finito-dimensional de \((H(\boldsymbol{s_1}), \ldots, H(\boldsymbol{s_n}))\) for a mesma de \((H(\boldsymbol{s_1} + \boldsymbol{u}), \ldots, H(\boldsymbol{s_n} + \boldsymbol{u}))\), para quaisquer \(\boldsymbol{s_i} \in \mathbb{G}\), \(i=1, \ldots, n\), e \(\boldsymbol{s_i} + \boldsymbol{u} \in \mathbb{G}\). Quando isto ocorre, tem-se invariância sob translações no espaço para todas as distribuições unidimensionais implicando em média e variância constantes, isto é, \[\begin{eqnarray} \mu(\boldsymbol{s}) = \mu \quad \mbox{e} \quad Var[H(\boldsymbol{s})] = V, \quad \forall \boldsymbol{s} \in \mathbb{G}. \end{eqnarray}\] Além disso, tem-se que a distribuição conjunta de \((H(\boldsymbol{s_i}), H(\boldsymbol{s_i} + \boldsymbol{u}))\), para quaisquer \(\boldsymbol{s_i} \in \mathbb{G}\) e \(\boldsymbol{s_i} + \boldsymbol{u} \in \mathbb{G}\), depende somente do vetor de separação \(\boldsymbol{u}\). Ao modelar um processo aleatório de forma paramétrica, é comum atribuir uma distribuição de probabilidade ao processo. Suponha que os parâmetros desta distribuição dependam da média e da variância do processo. Sendo assim, quando o processo é estritamente estacionário e o número de localizações é grande, há uma redução significativa no número de parâmetros que definem a distribuição do processo.

  • O processo é dito ser fracamente estacionário (ou estacionário de segunda ordem) se possuir as seguintes características:

    • a média do processo, \(\mu(\boldsymbol{s})\), não variar com as coordenadas \(\boldsymbol{s}\), ou seja, \[\mu(\boldsymbol{s}) = \mu, \quad \forall \boldsymbol{s} \in \mathbb{G},\]
    • a covariância entre duas localizações distintas depender apenas do vetor de separação \(\boldsymbol{u}\), ou seja, \[cov(H(\boldsymbol{s}), H(\boldsymbol{s}+\boldsymbol{u})) = C(\boldsymbol{u}), \quad \forall \boldsymbol{s}, \boldsymbol{s} + \boldsymbol{u} \in \mathbb{G}.\]
    • o valor esperado \(E[Y^2(\boldsymbol{s})] < \infty\), \(\forall \boldsymbol{s} \in \mathbb{G}\).

A função \(C(\cdot)\) é chamada de covariograma.

  • Um processo estritamente estacionário só é fracamente estacionário se a condição (iii) for atendida. Um processo fracamente estacionário costuma ser chamado apenas de processo estacionário.

  • Diz-se que o processo é intrinsecamente estacionário se possuir as seguintes características:

    • \(E[H(\boldsymbol{s} + \boldsymbol{u})] = E[H(\boldsymbol{s})]\)
    • \(Var[H(\boldsymbol{s}+\boldsymbol{u}) - H(\boldsymbol{s})] = 2\gamma(\boldsymbol{u})\), para todo \(\boldsymbol{s}, \boldsymbol{s} + \boldsymbol{u} \in \mathbb{G}\).

  • Chama-se a função \(2\gamma(\boldsymbol{u})\) de variograma e a função \(\gamma(\boldsymbol{u})\) de semivariograma.

  • O processo é denominado isotrópico, quando o variograma depender do vetor de separação \(\boldsymbol{u}\) apenas através do seu comprimento \(||\boldsymbol{u}||\). Caso contrário, é chamado de anisotrópico. Quando o processo é intrinsecamente estacionário e isotrópico, diz-se que o processo é homogêneo. Quando o processo não tem pelo menos uma destas propriedades, diz-se que o processo é heterogêneo.

  • Quando o processo é homogêneo, tem-se que a função de covariância entre \(H(\boldsymbol{s}+\boldsymbol{u})\) e \(H(\boldsymbol{s})\) depende apenas do comprimento do vetor \(\boldsymbol{u}\) e, portanto, a variância do processo é constante ao longo da região \(\mathbb{G}\), i.e.,\[Var[H(\boldsymbol{s})] = V, \quad \forall \boldsymbol{s} \in \mathbb{G}.\] Desta forma, podemos decompor a função de covariância em um produto de dois fatores, sendo um deles a variância e o outro uma função de correlação, i.e., \[\begin{eqnarray*} cov(H(\boldsymbol{s}), H(\boldsymbol{s}+\boldsymbol{u})) = C(\boldsymbol{u}) = V \rho(||\boldsymbol{u}||; \boldsymbol{\psi}), \end{eqnarray*}\] onde \(\rho(\cdot; \boldsymbol{\psi})\) é uma função de correlação válida em \(\mathbb{R}^p\) e \(\boldsymbol{\psi}\) é o conjunto de parâmetros desta função. Sendo assim, o variograma pode ser reescrito como \[\begin{eqnarray} 2\gamma(||\boldsymbol{u}||) &=& Var[H(\boldsymbol{s}+\boldsymbol{u})] + Var[H(\boldsymbol{s})] - 2cov(H(\boldsymbol{s}+\boldsymbol{u}), H(\boldsymbol{s})) \nonumber \\ &=& 2V\left[ 1 - \rho(||\boldsymbol{u}||; \boldsymbol{\psi})\right]. \end{eqnarray}\]

  • Ao representar graficamente um variograma de um processo homogêneo, há três medidas que podem ser visualizadas por ele: o patamar, o efeito pepita e o alcance. Por definição, quando o comprimento do vetor de separação é nulo, temos que o variograma é nulo, ou seja, \(2\gamma(0) = 0\). Além disso, espera-se que o processo aleatório esteja mais correlacionado à medida que a distância entre duas localizações tenda a zero. Sendo assim, espera-se que o variograma também tenda a zero para distâncias pequenas. Mas, na prática, isto dificilmente ocorre. Por isso, a medida obtida para o variograma quando a distância for nula, é chamada de efeito pepita e é considerada como um erro de medida. O alcance representa a distância mínima na qual a correlação entre duas localizações se anula. Quando a correlação se anula, tem-se que o variograma mantém-se constante. O patamar é o valor da variância quando a distância tende a infinito. A Figura abaixo apresenta um exemplo de um variograma com efeito pepita igual a \(0,4\), alcance igual a \(0,5\) e com patamar igual a \(1,6\).

  • Uma função de covariância só é considerada válida se for positiva semidefinida, ou seja, se para qualquer coleção finita de localizações \(\{\boldsymbol{s_1}, \ldots, \boldsymbol{s_n}\}\) e para quaisquer números reais \(\{a_1, \ldots, a_n\}\), então \[\begin{eqnarray} Var\left[ \sum_{i=1}^n{a_i (H(\boldsymbol{s_i})}\right] = \sum_{i=1}^n{\sum_{l=1}^n{a_i a_l Cov(H(\boldsymbol{s_i}), (H(\boldsymbol{s_l})) }} \geq 0. \end{eqnarray}\] Verificar se uma determinada função de covariância é válida é uma tarefa árdua. Por isso, muitas vezes, utiliza-se uma função de covariância conhecida na literatura. A tabela abaixo apresenta algumas funções de covariância válidas para processos homogêneos.

Exemplos de funções de covariância válidas para processos homogêneos onde \(u = ||\boldsymbol{u}||\), \(\tau^2, V \geq 0\) e \(\phi>0\). O parâmetro \(V\) representa o efeito pepita.

Nome Função de covariância
Esférica \(C(u) = \left\{\begin{array}{ll} \tau^2 \left[ 1 - 1,5 \phi^{-1} u + 0,5 (\phi^{-1} u)^3 \right], & \mbox{ se } 0 \leq u \leq \phi, \\ 0, & \mbox{ se } u \geq \phi, \\ V + \tau^2, & \mbox{se } u=0. \\ \end{array}\right .\)
Exponencial Potência \(C(u) = \left\{\begin{array}{ll} V + \tau^2, & \mbox{ se } u = 0, \\ \tau^2 \exp\left\{-u^{k} \phi^{-1}\right\}, & \mbox{ se } u>0, \end{array}\right .\) onde \(k \in (0,2]\).
Exponencial Potência/Quadrática \(C(u) = \left\{\begin{array}{ll} V + \tau^2, & \mbox{ se } u = 0, \\ \tau^2 \exp\left\{-u^2 \phi^{-2}\right\}, & \mbox{ se } u>0 \end{array}\right .\)
Matérn \(C(u) = \left\{\begin{array}{ll} V + \tau^2, & \mbox{ se } u = 0, \\ \frac{\tau^2}{2^{\lambda-1}\phi(\lambda)} \left( \frac{2\sqrt{\lambda} u}{\phi}\right) K_\lambda\left( \frac{2\sqrt{\lambda} u}{\phi}\right), & \mbox{ se } u>0 \end{array}\right .\) onde \(K_\lambda\) é a função modificada de Bessel do terceiro tipo de ordem \(\lambda\) e \(\Gamma(\cdot)\) é a função gama.
  • A função de covariância esférica é válida somente para dimensões \(p \leq 3\).
  • A função exponencial potência é válida em todas as dimensões \(p\) e é chamada de exponencial potência quadrática ou gaussiana, quando \(k = 2\), e de exponencial, quando \(k=1\). A função exponencial potência quadrática é utilizada para processos muito suaves. Note que a função exponencial potência produz somente correlações positivas.
  • A função Matérn tem como casos especiais a função exponencial, quando \(\lambda = 0,5\), e a exponencial potência quadrática, quando \(\lambda \rightarrow \infty\). Além disso, temos que o parâmetro \(\lambda\) controla a suavidade do processo.
  • Na função de covariância exponencial, a função de correlação, \(\exp\left\{-u \phi^{-1}\right\}\), nunca se anula. Por isso, o patamar só é alcançado assintoticamente e o alcance é infinito. Em casos como este, alguns autores (Banerjee, Gelfand, e Carlin (2015)) utilizam a noção de um alcance efetivo, que é a distância mínima que torna a correlação entre duas localizações distintas muito pequena. Em geral, considera-se que a correlação é muito pequena quando assume valores menores que \(0,05\).

Mais informações podem ser obtidas em Cressie (1993), Banerjee, Gelfand, e Carlin (2015) e Blangiardo e Cameletti (2015).

2 Modelo usando a função exponencial

\[\begin{eqnarray} Y(\boldsymbol{s}) &=& \boldsymbol{X}(\boldsymbol{s}) \boldsymbol{\beta} + H(\boldsymbol{s}) + e(\boldsymbol{s}), \end{eqnarray}\] onde \(\boldsymbol{X}(\boldsymbol{s})\) é um vetor de ordem \(K\) contendo variáveis explicativas, \(\boldsymbol{\beta}\) é um vetor coluna de ordem \(K\) representando os efeitos do intercepto e das covariáveis, \(H(\boldsymbol{s})\) é um efeito espacial e independente de \(e(\boldsymbol{s})\), que é um efeito aleatório.

Considere que \[\begin{eqnarray} e(\boldsymbol{s}) &\sim& N(0,V)\\ \boldsymbol{H} = (H(\boldsymbol{s}_1), \ldots, H(\boldsymbol{s}_n))^{'} &\sim& N_n(\boldsymbol{0}, \boldsymbol{\Sigma}) \end{eqnarray}\] sendo \(V \in \mathbb{R}^{+}\) e \(\boldsymbol{\Sigma}\) uma matriz de covariância na qual cada elemento é igual a \(\Sigma_{ij} = \tau^2 \exp\left(-\frac{d_{ij}}{\phi}\right)\), com \(\phi>0\).

Sendo assim, temos que o vetor de parâmetros desconhecidos é \[(\boldsymbol{\beta}, V, \tau^2, \phi).\]

Considere a seguinte distribuição a priori: \[\begin{eqnarray} p(\boldsymbol{\beta}, V, \tau^2, \phi) &=& p(\boldsymbol{\beta}) p(V) p(\tau^2) p(\phi)\\ \boldsymbol{\beta} &\sim& N_K(\boldsymbol{m}_\beta, \boldsymbol{V}_\beta)\\ V &\sim& GI(a_V, b_V)\\ \tau^2 &\sim& GI(a_\tau, b_\tau)\\ \phi &\sim& G(a_\phi, b_\phi) \end{eqnarray}\]

2.1 Simulação do conjunto de dados

Exercício: Simule um conjunto de dados do modelo acima. Para isso, siga os passos abaixo:

  • Considere que há 100 localizações.

  • Crie uma grade regular no quadrado unitário e para isso considere que a latitude e a longitude assumam valores entre 0 e 1.

  • Assuma que o intercepto vale 1 e que o efeito da covariável vale 5.

  • Gere uma covariável da normal padrão.

  • Assuma \(\phi = 1/6\).

  • Assuma \(V = 0,1\) e \(\tau^2=2\).

if(!require(MASS)){ install.packages("MASS"); require(MASS)}
if(!require(gstat)){ install.packages("gstat"); require(gstat)}
if(!require(sp)){ install.packages("sp"); require(sp)}
if(!require(fields)){ install.packages("fields"); require(fields)}
if(!require(spBayes)){ install.packages("spBayes"); require(spBayes)}

#Gerando pontos considerando uma grade regular no quadrado unitario
n   = 100 #total de regiões
lat = lon = NULL
lat = runif(n, 0, 1)
lon = runif(n, 0, 1)

#visualizando os pontos gerados
plot(lat, lon, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     xlab="latitude", ylab="longitude")

#calculando as distancias euclideanas entre as regioes
d = matrix(NA, n, n)
for(i in 1:n){
  for(j in 1:n){
    d[i,j] = sqrt( (lat[i] - lat[j])^2 + (lon[i]-lon[j])^2)
  }
}
#visualizando as distancias obtidas
hist(d[upper.tri(d)], freq=F, main="Histograma",
     xlab="distancias", ylab="densidade")

#Simulando os dados
beta  = c(1,5)
K     = length(beta)
x1    = rnorm(n)
x     = matrix(NA, n, K)
x[,1] = rep(1,n)
x[,2] = x1
V     = 0.1
e     = rnorm(n, 0, sqrt(V))
tau2  = 2 
phi   = 0.5/3#max(d)/6
Sigma = tau2 * exp(-d/phi)
H     = mvrnorm(1, rep(0,n), Sigma)
y     = x %*% beta + H + e
dados = data.frame(latitude=lat, longitude=lon, y=y, x=x1)

par(mfrow=c(2,2))
hist(y, freq=F, bty="n", cex.lab=1.4, cex.axis=1.4, lwd=2, main="Histograma",
     xlab="y", ylab="densidade")
hist(x %*% beta, freq=F, bty="n", cex.lab=1.4, cex.axis=1.4, lwd=2, main="Histograma",
     xlab="x %*% beta", ylab="densidade")
hist(H, freq=F, bty="n", cex.lab=1.4, cex.axis=1.4, lwd=2, main="Histograma",
     xlab="H", ylab="densidade")
hist(e, freq=F, bty="n", cex.lab=1.4, cex.axis=1.4, lwd=2, main="Histograma",
     xlab="e", ylab="densidade")

par(mfrow=c(1,2))
#Plotando o variograma
u     = seq(0, max(d), l=100)
vari  = function(u){return(2*tau2*(1-exp(-u/phi)))}
alcance_efetivo = -phi*log(0.05)
plot(u, vari(u), type="l", lwd=2, 
     bty="n", cex.lab=1.4, cex.axis=1.4,
     xlab="distancia", ylab="variograma")
abline(v=alcance_efetivo, lty=3, lwd=2)
abline(h=vari(alcance_efetivo), lty=3, lwd=2)
abline(h=2*tau2, lty=3, lwd=2, col="red")

#Plotando a funcao de correlacao
plot(u, exp(-u/phi), type="l", lwd=2, 
     bty="n", cex.lab=1.4, cex.axis=1.4,
     xlab="distancia", ylab="correlação")
abline(v=alcance_efetivo, lty=3, lwd=2)
abline(h=exp(-alcance_efetivo/phi), lty=3, lwd=2)

#visualizando a variavel observada
par(mfrow=c(1,1))
plot(lat, lon, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     xlab="latitude", ylab="longitude", cex=y)

##Ajustando um modelo linear normal bayesiano independente

#Ajustando um modelo linear normal bayesiano independente
#usando um pacote rstanarm
if(!require(rstanarm)){ install.packages("rstanarm"); require(rstanarm)}

modelo  = stan_glm(y ~ x -1, family = gaussian, verbose = FALSE)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.068 seconds (Warm-up)
## Chain 1:                0.048 seconds (Sampling)
## Chain 1:                0.116 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.041 seconds (Warm-up)
## Chain 2:                0.051 seconds (Sampling)
## Chain 2:                0.092 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.041 seconds (Warm-up)
## Chain 3:                0.049 seconds (Sampling)
## Chain 3:                0.09 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.043 seconds (Warm-up)
## Chain 4:                0.049 seconds (Sampling)
## Chain 4:                0.092 seconds (Total)
## Chain 4:
#comparando os valores verdadeiros com os estimados
#usando a media a posteriori
cbind(modelo$coefficients, beta)
##             beta
## x1 1.335670    1
## x2 5.137295    5
summary(modelo)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      y ~ x - 1
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 100
##  predictors:   2
## 
## Estimates:
##         mean   sd   10%   50%   90%
## x1    1.3    0.1  1.2   1.3   1.5  
## x2    5.1    0.1  5.0   5.1   5.3  
## sigma 1.3    0.1  1.2   1.3   1.5  
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 1.1    0.2  0.9   1.1   1.4  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## x1            0.0  1.0  3722 
## x2            0.0  1.0  3534 
## sigma         0.0  1.0  3479 
## mean_PPD      0.0  1.0  4035 
## log-posterior 0.0  1.0  1937 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
if(!require(bayesplot)){ install.packages("bayesplot"); require(bayesplot)}
#analisando a convergencia
mcmc_trace(modelo)

mcmc_acf(modelo)

#analisando os residuos
par(mfrow=c(1,2))
hist(modelo$residuals, freq=F, main="Histograma", xlab="residuos",
    ylab="densidade")
qqnorm(modelo$residuals, bty="n", lwd=2)
qqline(modelo$residuals, lwd=2, col="red")

2.2 Ajustando um modelo linear espacial usando função de covariância exponencial

# Ajustando um modelo linear espacial usando função de covariância exponencial

coordenadas = cbind(lat, lon)

starting <- list("phi"=3/0.5, "sigma.sq"=10, "tau.sq"=1)

tuning <- list("phi"=1/0.1, "sigma.sq"=0.1, "tau.sq"=0.1)

priors.1 <- list("beta.Norm"=list(rep(0,K), diag(1000,K)),
                 "phi.Unif"=c(3/1, 3/0.1), 
                 "sigma.sq.IG"=c(2, 2),
                 "tau.sq.IG"=c(2, 0.1))

priors.2 <- list("beta.Flat", "phi.Unif"=c(3/1, 3/0.1),
                 "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1))

cov.model <- "exponential"
n.samples <- 2000

m.1 <- spLM(y~x-1, coords=coordenadas, starting=starting,
            tuning=tuning, priors=priors.1, cov.model=cov.model,
            n.samples=n.samples, verbose=TRUE)
## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 100 observations.
## 
## Number of covariates 2 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Number of MCMC samples 2000.
## 
## Priors and hyperpriors:
##  beta normal:
##  mu: 0.000   0.000   
##  cov:
##  1000.000    0.000   
##  0.000   1000.000    
## 
##  sigma.sq IG hyperpriors shape=2.00000 and scale=2.00000
##  tau.sq IG hyperpriors shape=2.00000 and scale=0.10000
##  phi Unif hyperpriors a=3.00000 and b=30.00000
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Sampled: 100 of 2000, 5.00%
## Report interval Metrop. Acceptance rate: 21.00%
## Overall Metrop. Acceptance rate: 21.00%
## -------------------------------------------------
## Sampled: 200 of 2000, 10.00%
## Report interval Metrop. Acceptance rate: 11.00%
## Overall Metrop. Acceptance rate: 16.00%
## -------------------------------------------------
## Sampled: 300 of 2000, 15.00%
## Report interval Metrop. Acceptance rate: 9.00%
## Overall Metrop. Acceptance rate: 13.67%
## -------------------------------------------------
## Sampled: 400 of 2000, 20.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 11.75%
## -------------------------------------------------
## Sampled: 500 of 2000, 25.00%
## Report interval Metrop. Acceptance rate: 12.00%
## Overall Metrop. Acceptance rate: 11.80%
## -------------------------------------------------
## Sampled: 600 of 2000, 30.00%
## Report interval Metrop. Acceptance rate: 12.00%
## Overall Metrop. Acceptance rate: 11.83%
## -------------------------------------------------
## Sampled: 700 of 2000, 35.00%
## Report interval Metrop. Acceptance rate: 9.00%
## Overall Metrop. Acceptance rate: 11.43%
## -------------------------------------------------
## Sampled: 800 of 2000, 40.00%
## Report interval Metrop. Acceptance rate: 11.00%
## Overall Metrop. Acceptance rate: 11.38%
## -------------------------------------------------
## Sampled: 900 of 2000, 45.00%
## Report interval Metrop. Acceptance rate: 15.00%
## Overall Metrop. Acceptance rate: 11.78%
## -------------------------------------------------
## Sampled: 1000 of 2000, 50.00%
## Report interval Metrop. Acceptance rate: 16.00%
## Overall Metrop. Acceptance rate: 12.20%
## -------------------------------------------------
## Sampled: 1100 of 2000, 55.00%
## Report interval Metrop. Acceptance rate: 9.00%
## Overall Metrop. Acceptance rate: 11.91%
## -------------------------------------------------
## Sampled: 1200 of 2000, 60.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 11.75%
## -------------------------------------------------
## Sampled: 1300 of 2000, 65.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 11.31%
## -------------------------------------------------
## Sampled: 1400 of 2000, 70.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 11.07%
## -------------------------------------------------
## Sampled: 1500 of 2000, 75.00%
## Report interval Metrop. Acceptance rate: 5.00%
## Overall Metrop. Acceptance rate: 10.67%
## -------------------------------------------------
## Sampled: 1600 of 2000, 80.00%
## Report interval Metrop. Acceptance rate: 3.00%
## Overall Metrop. Acceptance rate: 10.19%
## -------------------------------------------------
## Sampled: 1700 of 2000, 85.00%
## Report interval Metrop. Acceptance rate: 9.00%
## Overall Metrop. Acceptance rate: 10.12%
## -------------------------------------------------
## Sampled: 1800 of 2000, 90.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 10.33%
## -------------------------------------------------
## Sampled: 1900 of 2000, 95.00%
## Report interval Metrop. Acceptance rate: 15.00%
## Overall Metrop. Acceptance rate: 10.58%
## -------------------------------------------------
## Sampled: 2000 of 2000, 100.00%
## Report interval Metrop. Acceptance rate: 5.00%
## Overall Metrop. Acceptance rate: 10.30%
## -------------------------------------------------
m.2 <- spLM(y~x-1, coords=coordenadas, starting=starting,
            tuning=tuning, priors=priors.2, cov.model=cov.model,
            n.samples=n.samples, verbose=TRUE)
## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 100 observations.
## 
## Number of covariates 2 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Number of MCMC samples 2000.
## 
## Priors and hyperpriors:
##  beta flat.
##  sigma.sq IG hyperpriors shape=2.00000 and scale=2.00000
##  tau.sq IG hyperpriors shape=2.00000 and scale=0.10000
##  phi Unif hyperpriors a=3.00000 and b=30.00000
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Sampled: 100 of 2000, 5.00%
## Report interval Metrop. Acceptance rate: 17.00%
## Overall Metrop. Acceptance rate: 17.00%
## -------------------------------------------------
## Sampled: 200 of 2000, 10.00%
## Report interval Metrop. Acceptance rate: 1.00%
## Overall Metrop. Acceptance rate: 9.00%
## -------------------------------------------------
## Sampled: 300 of 2000, 15.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 10.33%
## -------------------------------------------------
## Sampled: 400 of 2000, 20.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 9.25%
## -------------------------------------------------
## Sampled: 500 of 2000, 25.00%
## Report interval Metrop. Acceptance rate: 11.00%
## Overall Metrop. Acceptance rate: 9.60%
## -------------------------------------------------
## Sampled: 600 of 2000, 30.00%
## Report interval Metrop. Acceptance rate: 20.00%
## Overall Metrop. Acceptance rate: 11.33%
## -------------------------------------------------
## Sampled: 700 of 2000, 35.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 10.57%
## -------------------------------------------------
## Sampled: 800 of 2000, 40.00%
## Report interval Metrop. Acceptance rate: 15.00%
## Overall Metrop. Acceptance rate: 11.12%
## -------------------------------------------------
## Sampled: 900 of 2000, 45.00%
## Report interval Metrop. Acceptance rate: 3.00%
## Overall Metrop. Acceptance rate: 10.22%
## -------------------------------------------------
## Sampled: 1000 of 2000, 50.00%
## Report interval Metrop. Acceptance rate: 12.00%
## Overall Metrop. Acceptance rate: 10.40%
## -------------------------------------------------
## Sampled: 1100 of 2000, 55.00%
## Report interval Metrop. Acceptance rate: 7.00%
## Overall Metrop. Acceptance rate: 10.09%
## -------------------------------------------------
## Sampled: 1200 of 2000, 60.00%
## Report interval Metrop. Acceptance rate: 4.00%
## Overall Metrop. Acceptance rate: 9.58%
## -------------------------------------------------
## Sampled: 1300 of 2000, 65.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 9.62%
## -------------------------------------------------
## Sampled: 1400 of 2000, 70.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 9.36%
## -------------------------------------------------
## Sampled: 1500 of 2000, 75.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 9.67%
## -------------------------------------------------
## Sampled: 1600 of 2000, 80.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 9.94%
## -------------------------------------------------
## Sampled: 1700 of 2000, 85.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 9.82%
## -------------------------------------------------
## Sampled: 1800 of 2000, 90.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 9.72%
## -------------------------------------------------
## Sampled: 1900 of 2000, 95.00%
## Report interval Metrop. Acceptance rate: 7.00%
## Overall Metrop. Acceptance rate: 9.58%
## -------------------------------------------------
## Sampled: 2000 of 2000, 100.00%
## Report interval Metrop. Acceptance rate: 12.00%
## Overall Metrop. Acceptance rate: 9.70%
## -------------------------------------------------
burn.in <- 100

##recover beta and spatial random effects
m.1 <- spRecover(m.1, start=burn.in, verbose=FALSE)
m.2 <- spRecover(m.2, start=burn.in, verbose=FALSE)

cbind(c(tau2, V, 1/phi),
      round(summary(m.1$p.theta.recover.samples)$quantiles[,c(3,1,5)],2),
      round(summary(m.2$p.theta.recover.samples)$quantiles[,c(3,1,5)],2))
##                50% 2.5% 97.5%   50% 2.5% 97.5%
## sigma.sq 2.0  1.55 1.08  2.27  1.56 1.01  2.34
## tau.sq   0.1  0.05 0.02  0.21  0.07 0.02  0.17
## phi      6.0 10.42 5.50 17.02 10.31 5.62 16.72
cbind(beta,
      round(summary(m.1$p.beta.recover.samples)$quantiles[,c(3,1,5)],2),
      round(summary(m.2$p.beta.recover.samples)$quantiles[,c(3,1,5)],2))
##    beta  50% 2.5% 97.5%  50% 2.5% 97.5%
## x1    1 1.25 0.67  1.87 1.25 0.71  1.87
## x2    5 5.06 4.87  5.24 5.05 4.87  5.24
library(coda)
m.1.w.summary <- summary(mcmc(t(m.1$p.w.recover.samples)))$quantiles[,c(3,1,5)]
m.2.w.summary <- summary(mcmc(t(m.2$p.w.recover.samples)))$quantiles[,c(3,1,5)]

par(mfrow=c(1,1))
plot(H, m.1.w.summary[,1], xlab="observado", ylab="ajustado",
     xlim=range(H), ylim=range(m.1.w.summary), main="Efeito espacial H")
arrows(H, m.1.w.summary[,1], H, m.1.w.summary[,2], length=0.02, angle=90)
arrows(H, m.1.w.summary[,1], H, m.1.w.summary[,3], length=0.02, angle=90)
lines(range(H), range(H))

points(H, m.2.w.summary[,1], col="blue", pch=19, cex=0.5)
arrows(H, m.2.w.summary[,1], H, col="blue", m.2.w.summary[,2], length=0.02, angle=90)
arrows(H, m.2.w.summary[,1], H, col="blue", m.2.w.summary[,3], length=0.02, angle=90)

###########################
##Predictive process model
###########################
m.1 <- spLM(y~x-1, coords=coordenadas, knots=c(6,6,0.1), starting=starting,
            tuning=tuning, priors=priors.1, cov.model=cov.model,
            n.samples=n.samples, verbose=TRUE)
## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 100 observations.
## 
## Number of covariates 2 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Using modified predictive process with 36 knots.
## 
## Number of MCMC samples 2000.
## 
## Priors and hyperpriors:
##  beta normal:
##  mu: 0.000   0.000   
##  cov:
##  1000.000    0.000   
##  0.000   1000.000    
## 
##  sigma.sq IG hyperpriors shape=2.00000 and scale=2.00000
##  tau.sq IG hyperpriors shape=2.00000 and scale=0.10000
##  phi Unif hyperpriors a=3.00000 and b=30.00000
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Sampled: 100 of 2000, 5.00%
## Report interval Metrop. Acceptance rate: 27.00%
## Overall Metrop. Acceptance rate: 27.00%
## -------------------------------------------------
## Sampled: 200 of 2000, 10.00%
## Report interval Metrop. Acceptance rate: 15.00%
## Overall Metrop. Acceptance rate: 21.00%
## -------------------------------------------------
## Sampled: 300 of 2000, 15.00%
## Report interval Metrop. Acceptance rate: 20.00%
## Overall Metrop. Acceptance rate: 20.67%
## -------------------------------------------------
## Sampled: 400 of 2000, 20.00%
## Report interval Metrop. Acceptance rate: 16.00%
## Overall Metrop. Acceptance rate: 19.50%
## -------------------------------------------------
## Sampled: 500 of 2000, 25.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 17.60%
## -------------------------------------------------
## Sampled: 600 of 2000, 30.00%
## Report interval Metrop. Acceptance rate: 19.00%
## Overall Metrop. Acceptance rate: 17.83%
## -------------------------------------------------
## Sampled: 700 of 2000, 35.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 16.71%
## -------------------------------------------------
## Sampled: 800 of 2000, 40.00%
## Report interval Metrop. Acceptance rate: 6.00%
## Overall Metrop. Acceptance rate: 15.38%
## -------------------------------------------------
## Sampled: 900 of 2000, 45.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 15.22%
## -------------------------------------------------
## Sampled: 1000 of 2000, 50.00%
## Report interval Metrop. Acceptance rate: 21.00%
## Overall Metrop. Acceptance rate: 15.80%
## -------------------------------------------------
## Sampled: 1100 of 2000, 55.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 15.55%
## -------------------------------------------------
## Sampled: 1200 of 2000, 60.00%
## Report interval Metrop. Acceptance rate: 16.00%
## Overall Metrop. Acceptance rate: 15.58%
## -------------------------------------------------
## Sampled: 1300 of 2000, 65.00%
## Report interval Metrop. Acceptance rate: 18.00%
## Overall Metrop. Acceptance rate: 15.77%
## -------------------------------------------------
## Sampled: 1400 of 2000, 70.00%
## Report interval Metrop. Acceptance rate: 19.00%
## Overall Metrop. Acceptance rate: 16.00%
## -------------------------------------------------
## Sampled: 1500 of 2000, 75.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 15.80%
## -------------------------------------------------
## Sampled: 1600 of 2000, 80.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 15.44%
## -------------------------------------------------
## Sampled: 1700 of 2000, 85.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 15.35%
## -------------------------------------------------
## Sampled: 1800 of 2000, 90.00%
## Report interval Metrop. Acceptance rate: 11.00%
## Overall Metrop. Acceptance rate: 15.11%
## -------------------------------------------------
## Sampled: 1900 of 2000, 95.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 14.74%
## -------------------------------------------------
## Sampled: 2000 of 2000, 100.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 14.65%
## -------------------------------------------------
m.2 <- spLM(y~x-1, coords=coordenadas, knots=c(6,6,0.1), starting=starting,
            tuning=tuning, priors=priors.2, cov.model=cov.model,
            n.samples=n.samples, verbose=TRUE)
## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 100 observations.
## 
## Number of covariates 2 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Using modified predictive process with 36 knots.
## 
## Number of MCMC samples 2000.
## 
## Priors and hyperpriors:
##  beta flat.
##  sigma.sq IG hyperpriors shape=2.00000 and scale=2.00000
##  tau.sq IG hyperpriors shape=2.00000 and scale=0.10000
##  phi Unif hyperpriors a=3.00000 and b=30.00000
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Sampled: 100 of 2000, 5.00%
## Report interval Metrop. Acceptance rate: 25.00%
## Overall Metrop. Acceptance rate: 25.00%
## -------------------------------------------------
## Sampled: 200 of 2000, 10.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 16.50%
## -------------------------------------------------
## Sampled: 300 of 2000, 15.00%
## Report interval Metrop. Acceptance rate: 10.00%
## Overall Metrop. Acceptance rate: 14.33%
## -------------------------------------------------
## Sampled: 400 of 2000, 20.00%
## Report interval Metrop. Acceptance rate: 15.00%
## Overall Metrop. Acceptance rate: 14.50%
## -------------------------------------------------
## Sampled: 500 of 2000, 25.00%
## Report interval Metrop. Acceptance rate: 22.00%
## Overall Metrop. Acceptance rate: 16.00%
## -------------------------------------------------
## Sampled: 600 of 2000, 30.00%
## Report interval Metrop. Acceptance rate: 18.00%
## Overall Metrop. Acceptance rate: 16.33%
## -------------------------------------------------
## Sampled: 700 of 2000, 35.00%
## Report interval Metrop. Acceptance rate: 18.00%
## Overall Metrop. Acceptance rate: 16.57%
## -------------------------------------------------
## Sampled: 800 of 2000, 40.00%
## Report interval Metrop. Acceptance rate: 16.00%
## Overall Metrop. Acceptance rate: 16.50%
## -------------------------------------------------
## Sampled: 900 of 2000, 45.00%
## Report interval Metrop. Acceptance rate: 17.00%
## Overall Metrop. Acceptance rate: 16.56%
## -------------------------------------------------
## Sampled: 1000 of 2000, 50.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 16.30%
## -------------------------------------------------
## Sampled: 1100 of 2000, 55.00%
## Report interval Metrop. Acceptance rate: 14.00%
## Overall Metrop. Acceptance rate: 16.09%
## -------------------------------------------------
## Sampled: 1200 of 2000, 60.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 15.83%
## -------------------------------------------------
## Sampled: 1300 of 2000, 65.00%
## Report interval Metrop. Acceptance rate: 13.00%
## Overall Metrop. Acceptance rate: 15.62%
## -------------------------------------------------
## Sampled: 1400 of 2000, 70.00%
## Report interval Metrop. Acceptance rate: 8.00%
## Overall Metrop. Acceptance rate: 15.07%
## -------------------------------------------------
## Sampled: 1500 of 2000, 75.00%
## Report interval Metrop. Acceptance rate: 16.00%
## Overall Metrop. Acceptance rate: 15.13%
## -------------------------------------------------
## Sampled: 1600 of 2000, 80.00%
## Report interval Metrop. Acceptance rate: 26.00%
## Overall Metrop. Acceptance rate: 15.81%
## -------------------------------------------------
## Sampled: 1700 of 2000, 85.00%
## Report interval Metrop. Acceptance rate: 21.00%
## Overall Metrop. Acceptance rate: 16.12%
## -------------------------------------------------
## Sampled: 1800 of 2000, 90.00%
## Report interval Metrop. Acceptance rate: 18.00%
## Overall Metrop. Acceptance rate: 16.22%
## -------------------------------------------------
## Sampled: 1900 of 2000, 95.00%
## Report interval Metrop. Acceptance rate: 24.00%
## Overall Metrop. Acceptance rate: 16.63%
## -------------------------------------------------
## Sampled: 2000 of 2000, 100.00%
## Report interval Metrop. Acceptance rate: 17.00%
## Overall Metrop. Acceptance rate: 16.65%
## -------------------------------------------------
burn.in <- 0.5*n.samples

round(summary(window(m.1$p.beta.samples, start=burn.in))$quantiles[,c(3,1,5)],2)
##     50% 2.5% 97.5%
## x1 1.39 0.63  2.35
## x2 5.11 4.91  5.31
round(summary(window(m.2$p.beta.samples, start=burn.in))$quantiles[,c(3,1,5)],2)
##     50% 2.5% 97.5%
## x1 1.43 0.56  2.46
## x2 5.11 4.89  5.31
round(summary(window(m.1$p.theta.samples, start=burn.in))$quantiles[,c(3,1,5)],2)
##           50% 2.5% 97.5%
## sigma.sq 1.75 1.17  2.73
## tau.sq   0.07 0.03  0.18
## phi      5.48 3.27 10.27
round(summary(window(m.2$p.theta.samples, start=burn.in))$quantiles[,c(3,1,5)],2)
##           50% 2.5% 97.5%
## sigma.sq 1.88 1.26  2.96
## tau.sq   0.04 0.01  0.28
## phi      5.07 3.10  8.56

Referências

Banerjee, Sudipto, Alan E. Gelfand, e Bradley P. Carlin. 2015. Hierarchical Modeling and Analysis for Spatial Data. Chapman & Hall/CRC.
Blangiardo, Marta, e Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.
Cressie, Noel A. C. 1993. Statistics for Spatial Data. John Wiley & Sons.