Modelagem Espacial Bayesiana no R
Patrícia Lusié Velozo da Costa
27, 29 e 31 de março de 2023
1ª Escola de Verão - GET/IME/UFF
1 Estatística Espacial
É usada em muitas áreas, incluindo geografia, ciências ambientais, geologia, biologia, economia, demografia, criminologia, epidemiologia, engenharia, planejamento urbano e muitas outras.
Permite entender melhor como as variáveis se comportam em diferentes locais e como elas interagem umas com as outras em diferentes escalas espaciais. Isso pode ser usado para fazer previsões em áreas onde os dados são limitados, melhorar o planejamento e tomada de decisão em políticas públicas, bem como contribuir para o desenvolvimento de teorias em diferentes áreas do conhecimento.
Referências: Cressie (1993), Banerjee, Gelfand, e Carlin (2015) e Blangiardo e Cameletti (2015).
1.1 Tipos de dados espaciais
Costuma-se classificar em 3 tipos de dados espaciais:
- Dados de Área.
Nos dados de área, a variável de interesse em cada unidade elementar é um valor agregado em uma unidade de área com limites bem definidos como, por exemplo, bairros e cidades.
Dentro de cada área, supõe-se haver homogeneidade interna, ou seja, mudanças importantes só ocorrem nos limites.
- Superfícies Contínuas ou Geoestatística.
- São realizações de um processo espacial caracterizado por um índice espacial \(s\) que varia continuamente no domínio fixo.
- Eventos ou Padrões Pontuais.
- A incerteza está na localização de ocorrência de um determinado fenômeno como, por exemplo, a localização de crimes e a ocorrências de determinada doença.
Casos da doença pé e boca na região norte da Cumbria em 2001.
Fonte: Blangiardo e Cameletti (2015).
Na aula de hoje nos concentraremos em modelos para dados de área.
2 Modelo Gaussiano Espacial Condicional Autoregressivo (CAR) para dados de área
Considere o seguinte modelo: \[\begin{eqnarray} Y_i &=& \boldsymbol{x}_i^T \boldsymbol{\beta} + \phi_i + e_i, \quad e_i \sim N(0,V), \end{eqnarray}\] sendo
- \(\boldsymbol{x}_i\) o vetor de covariáveis da \(i\)-ésima região,
- \(\boldsymbol{\beta}\) o vetor de efeitos das covariáveis,
- \(V\) um parâmetro de variância e
- \(\phi_i\) o efeito aleatório espacial da \(i\)-ésima região.
Seja \(\boldsymbol{\phi}^T = (\phi_1, \ldots, \phi_n)\). A priori, considere que \((\boldsymbol{\beta}, V, \boldsymbol{\phi})\) são desconhecidos e que
\[\begin{eqnarray} p(\boldsymbol{\beta}, V, \boldsymbol{\phi}) &=& p(\boldsymbol{\beta})p(V)p(\boldsymbol{\phi}). \end{eqnarray}\]Para modelar o efeito espacial considere a distribuição a priori condicional autoregressiva (CAR), proposta por Leroux, Lei, e Breslow (2000): \[\begin{eqnarray} \phi_i | \boldsymbol{\phi}_{-i}, \boldsymbol{W}, \tau^2, \rho &\sim& N\left(\frac{\rho \sum_{j=1}^n{w_{ij} \phi_j}}{\rho \sum_{j=1}^n{w_{ij}} + 1 - \rho}, \frac{\tau^2}{\rho \sum_{j=1}^n{w_{ij}} + 1 - \rho}\right)\\ \tau^2 &\sim& GI(a_\tau, b_\tau)\\ \rho &\sim& Unif(0,1) \end{eqnarray}\] sendo
- \(\tau^2\) a variância do efeito espacial;
- \(\rho\) um parâmetro de correlação;
- \(\boldsymbol{W}\) uma matriz de vizinhança na qual cada elemento assume o valor 1 quando diferentes localizações são vizinhas e 0 caso contrário.
Quando \(\rho = 1\), tem-se o modelo CAR intrínseco.
Quando \(\rho = 0\), tem-se o modelo independente.
Tem-se que a distribuição conjunta dos efeitos espaciais é dada pela seguinte expressão \[\begin{eqnarray} \boldsymbol{\phi} \sim N_n\left(\boldsymbol{0}, \tau^2 \left(\rho(\boldsymbol{D} - \boldsymbol{W}) + (1-\rho) \boldsymbol{I} \right)^{-1}\right), \end{eqnarray}\] sendo \(\boldsymbol{I}\) a matriz identidade e \(\boldsymbol{D}\) uma matriz diagonal com o i-ésimo da diagonal igual ao número de vizinhos da região \(i\).
Pacote CARBayes
2.1 Ajustando um modelo independente
Quando \(\phi_i = 0, \forall i\), tem-se o modelo independente como caso particular do modelo apresentado acima.
Para aprendermos a usar o pacote, simularemos dados de um modelo de regressão linear e usaremos a função S.glm para ajustar este modelo.
Pode se incluir uma coleção de covariáveis espaciais neste modelo. Desta forma, o modelo não foi estruturado para ter dependência espacial mas dependendo das covariáveis pode capturar algum comportamento espacial.
Mais informações sobre modelos de regressão linear podem ser obtidas em Agresti (2015).
Usaremos os seguintes pacotes:
#Instalando (caso seja necessário e chamando pacotes necessários)
#Pacote para o ajuste do modelo espacial CAR
if(!require(CARBayes)){ install.packages("CARBayes"); require(CARBayes)}
#análise de convergência do MCMC
if(!require(coda)){ install.packages("coda"); require(coda)}
#Construção e customização de gráficos
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
#Permite o acesso aos shapefiles de diferentes unidades geográficas brasileiras.
if(!require(geobr)){ install.packages("geobr"); require(geobr)}
#gridExtra: colocar diferentes graficos numa mesma tela
if(!require(gridExtra)){ install.packages("gridExtra"); require(gridExtra)}
#manipulação e análise de dados espaciais
if(!require(sf)){ install.packages("sf"); require(sf)}
#classes e métodos para dados espaciais
if(!require(sp)){ install.packages("sp"); require(sp)}
#construção das matrizes de vizinhança
if(!require(spdep)){ install.packages("spdep"); require(spdep)}
#manipulação e análise de dados
if(!require(tidyverse)){ install.packages("tidyverse"); require(tidyverse)}- Simularemos dados considerando que as regiões de interesse sejam os municípios do Estado do RJ:
#definindo a regiao
Regiao = read_municipality(code_muni = 33,year=2020,
simplified = T, showProgress = F)
#pacote geobr
#Plotando as regiões de interesse
# Removendo eixos dos grafico
no_axis <- theme(axis.title=element_blank(), #pacote ggplot2
axis.text=element_blank(),
axis.ticks=element_blank())
ggplot() +
geom_sf(data=Regiao, color= "black") +
#geom_sf_text(data = Regiao, aes(label = name_muni),
# color="black",size=3)+
labs(subtitle="Municípios do RJ", size=8) +
theme_minimal() +
no_axis
Usaremos a seguinte matriz de vizinhança: \[\begin{eqnarray} W_{ij} &=&\left\{\begin{array}{ll} 0 & \mbox{ se $i=j$ e se as regiões $i$ e $j$ não tiverem qualquer tipo de fronteira}\\ 1 & \mbox{ se as regiões $i$ e $j$ tiverem qualquer tipo de fronteira} \end{array} \right . \end{eqnarray}\]
#Criando a matriz de vizinhanca usando o pacote Geobr
sf.uf <- sf:::as_Spatial(Regiao$geom) #Mudando de arquivo tipo SF para SP (pacote sf e sp)
W.nb <- spdep::poly2nb(sf.uf, queen=T) #Calculando os vizinhos pela contiguidade Queen
W.list <- nb2listw(W.nb, style="B")
W <- nb2mat(W.nb, style="B")- Simulando um conjunto de dados:
#####################################################
#Simulando um valor observado para cada uma das regiões
#sob o modelo de regressão linear
#####################################################
dados = Regiao
#Gerando um conjunto de dados
n = nrow(W)
V = 1
beta = c(1.5,2) #intercepto + efeito das covariaveis
x = matrix(NA,nrow=n,ncol=2) #covariaveis
x[,1] = rep(1, n)
x[,2] = rnorm(n, 0, 1)
dados$x = x[,2]
mu = x %*% beta
y = mu + rnorm(n,0,sqrt(V))
dados$y = y
#plotando a variavel resposta
hist(dados$y, freq=F, main="Histograma", cex.lab=1.4, cex.axis=1.4, ylab="densidade",
xlab="y")
#comparando a variável resposta e a covariável
#plotando a variável resposta
g1=ggplot() +
geom_sf(data = dados, aes(fill=y), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
#plotando a covariável
g2=ggplot() +
geom_sf(data = dados, aes(fill=x), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
gridExtra::grid.arrange(g1, g2, top = '',ncol=2)
#gráfico de dispersão
plot(dados$x, dados$y, cex.lab=1.4, cex.axis=1.4, ylab="y", xlab="x", bty="n",lwd=2)
lines(dados$x, mu, col="red", lwd=2)
#######################################################################
# Ajustando o Modelo independente usando o pacote CARBayes
#######################################################################
#Definindo as variaveis necessarias para o MCMC
n.aquec <- 100 #numero de iteracoes para o aquecimento
n.espac <- 1 #numero do espacamento
n.post <- 1000 #tamanho da amostra
n.ite <- n.aquec + n.post*n.espac #total de iteracoes
#Definindo os hiperparametros da distribuicao a priori
V_beta = 1000
a_V = 2
b_V = 1
K = 2 #intercepto+1 covariável
# Ajustando o modelo via CARBayes
ajuste_MRL_CAR = S.glm(formula=dados$y ~ dados$x,
family="gaussian",
burnin=n.aquec, n.sample=n.ite,
thin = n.espac, prior.mean.beta=rep(0,K),
prior.var.beta=rep(V_beta,K),
prior.nu2 = c(a_V, b_V), verbose=FALSE)
#Avaliando as estimativas pontuais e intervalares
#e comparando com os valores verdadeiros
cbind( c(beta, V), ajuste_MRL_CAR$summary.results)## Mean 2.5% 97.5% n.sample % accept n.effective Geweke.diag
## (Intercept) 1.5 1.5508 1.3359 1.7529 1000 100 1000 1.0
## dados$x 2.0 1.9851 1.7631 2.1945 1000 100 1000 0.1
## nu2 1.0 0.9905 0.7109 1.3170 1000 100 1000 0.7#Alocando as amostras a posteriori
V_post = ajuste_MRL_CAR$samples$nu2
beta_post = ajuste_MRL_CAR$samples$beta
#Analisando convergencia
par(mar=c(4,6,4,4),mfrow=c(2,3))
# Traços a posteriori
traceplot(beta_post[,1],ylab=expression(beta[1]),xlab="iterações",
cex.lab=1.8,cex.axis=1.6)
traceplot(beta_post[,2],ylab=expression(beta[2]),xlab="iterações",
cex.lab=1.8,cex.axis=1.6)
traceplot(V_post,xlab="iterações",ylab=expression(V),
cex.lab=1.8,cex.axis=1.6)
# Autocorrelação
acf(beta_post[,1],bty="n",ylab=expression(beta[1]),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(beta_post[,2],bty="n",ylab=expression(beta[2]),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(V_post,bty="n",ylab="V",main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
# Histograma
par(mfrow=c(1,3))
hist(beta_post[,1],freq=F,bty="n",ylab="densidade",xlab= expression(beta[1]),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=beta[1],lwd=2,col="red")
abline(v=mean(beta_post[,1]),lwd=2,col="blue")
abline(v=quantile(beta_post[,1],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(beta_post[,1],0.975),lwd=2,lty=3,col="blue")
hist(beta_post[,2],freq=F,bty="n",ylab="densidade",xlab= expression(beta[2]),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=beta[2],lwd=2,col="red")
abline(v=mean(beta_post[,2]),lwd=2,col="blue")
abline(v=quantile(beta_post[,2],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(beta_post[,2],0.975),lwd=2,lty=3,col="blue")
hist(V_post,freq=F,bty="n",ylab="densidade",xlab="V",
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=V,lwd=2,col="red")
abline(v=mean(V_post),lwd=2,col="blue")
abline(v=quantile(V_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(V_post,0.975),lwd=2,lty=3,col="blue")
#analisando os residuos do modelo
residuos_CAR = unname(ajuste_MRL_CAR$residuals$response)
#Histograma e qqplot dos resíduos do modelo
par(mar=c(6,5,6,2),mfrow=c(1,2))
hist(ajuste_MRL_CAR$residuals$response,freq=F,cex.lab=1.8,cex.axis=1.6,ylab="densidade",xlab="resíduos",main="")
qqnorm(ajuste_MRL_CAR$residuals$response,cex.lab=1.4,cex.axis=1.4,ylab="quantis amostrais",xlab="quantis teóricos",main="",lwd=2,bty="n")
qqline(ajuste_MRL_CAR$residuals$response,col="red",lwd=2)
#analisando normalidade dos residuos
shapiro.test(ajuste_MRL_CAR$residuals$response)##
## Shapiro-Wilk normality test
##
## data: ajuste_MRL_CAR$residuals$response
## W = 0.99213, p-value = 0.8654ks.test(ajuste_MRL_CAR$residuals$response, "pnorm",
mean(ajuste_MRL_CAR$residuals$response), sd(ajuste_MRL_CAR$residuals$response))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: ajuste_MRL_CAR$residuals$response
## D = 0.063456, p-value = 0.8296
## alternative hypothesis: two-sided#Calculando Indice de Moran para os resíduos
moran(ajuste_MRL_CAR$residuals$response, W.list, n = nrow(dados), S0 = Szero(W.list))$I## [1] 0.03544861#H0: nao ha dependencia espacial
#H1: ha dependencia espacial
IM=moran.test(ajuste_MRL_CAR$residuals$response, W.list)
IM##
## Moran I test under randomisation
##
## data: ajuste_MRL_CAR$residuals$response
## weights: W.list
##
## Moran I statistic standard deviate = 0.71592, p-value = 0.237
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.035448611 -0.010989011 0.004207375alfa = 0.05
pvalor = IM$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("Não rejeito a hipotese de independencia espacial.")
}## [1] "Não rejeito a hipotese de independencia espacial."2.2 Ajustando um modelo Gaussiano espacial
#####################################################
#Simulando um valor observado para cada uma das regiões
#sob o modelo gaussiano com distribuição a priori CAR Leroux
#####################################################
dados2 = Regiao
#Gerando um conjunto de dados
n = nrow(W)
beta = 2 #efeito da covariavel
K = length(beta)
x = rnorm(n, 0, sqrt(1))
I = diag(n) #matriz identidade
tau2 = 1 #variancia do termo espacial
Nviz = apply(W,1,sum) #total de vizinhos da regiao i
D = diag(Nviz)
rho = 0.9 #parametro de correlacao espacial
Sigma = tau2*solve(rho * (D-W) + (1-rho)*I) #matriz de covariancia do efeito espacial
phi = mvrnorm(1, rep(0,n), Sigma) #efeito espacial
mu = x * beta + phi
V = 1
y = mu + rnorm(1, 0, sqrt(V))
#Analisando a correlacao gerada para o efeito espacial phi
Corre = matrix(NA,n,n)
for(i in 1:n){
for(j in 1:n){
Corre[i,j] = Sigma[i,j] / sqrt(Sigma[i,i] * Sigma[j,j])
}
}
Corre[1,]## [1] 1.00000000 0.02420996 0.04913556 0.13427051 0.02201420 0.03527878
## [7] 0.33184532 0.34227453 0.20325271 0.04822478 0.02219614 0.03791270
## [13] 0.07929296 0.02843982 0.02980262 0.10845376 0.03048499 0.04800338
## [19] 0.02335836 0.06482203 0.05045869 0.03643402 0.04588393 0.06111380
## [25] 0.23054154 0.24508352 0.11226258 0.03312123 0.06686157 0.42544244
## [31] 0.02512363 0.03291333 0.02413638 0.14373842 0.25353661 0.01693958
## [37] 0.04598016 0.03884840 0.13841521 0.61412481 0.05024552 0.28880741
## [43] 0.20673266 0.26813075 0.02172962 0.01763805 0.18327435 0.03648050
## [49] 0.07411111 0.28608900 0.31706504 0.18969377 0.65087090 0.19324805
## [55] 0.19203020 0.28782698 0.38689201 0.01334828 0.22988668 0.24724755
## [61] 0.21718428 0.02684013 0.22898777 0.05937648 0.58605903 0.19634843
## [67] 0.03299066 0.25093637 0.03447170 0.02974885 0.01785759 0.03191569
## [73] 0.04358948 0.01785759 0.19661889 0.02314317 0.12364116 0.03527878
## [79] 0.03866168 0.09106240 0.04906960 0.27783655 0.05732230 0.08779118
## [85] 0.05702989 0.10929984 0.04905044 0.12796435 0.27110850 0.01560060
## [91] 0.27437716 0.36671324 aux = Corre
diag(aux) = NA
max(aux, na.rm=T)## [1] 0.7067318 min(aux, na.rm=T)## [1] 0.008688009 length(which(aux>0.5))## [1] 492 par(mfrow=c(1,1))
hist(aux,main="Histograma", xlab="correlação de phi", ylab="densidade",
cex.lab=1.4, cex.axis=1.4)
#plotando a variavel resposta
hist(y, freq=F, main="Histograma", cex.lab=1.4, cex.axis=1.4, ylab="densidade",
xlab="y")
dados2$y = y
dados2$x = x
#comparando a variável resposta e a covariável
#plotando a variável resposta
g1=ggplot() +
geom_sf(data = dados2, aes(fill=y), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
#plotando a covariável
g2=ggplot() +
geom_sf(data = dados2, aes(fill=x), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
gridExtra::grid.arrange(g1, g2, top = '',ncol=2)
#gráfico de dispersão
plot(dados2$x, dados2$y, cex.lab=1.4, cex.axis=1.4, ylab="y", xlab="x", bty="n",lwd=2)
#Analisando phi
hist(phi, freq=F, main="Histograma", xlab=expression(phi),
ylab="densidade", cex.lab=1.4, cex.axis=1.4)
#atraves do Indice de Moran
alfa=0.05
pvalor = moran.test(phi, W.list)$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("NC#o rejeito a hipotese de independencia espacial.")
}## [1] "Rejeito a hipotese de independencia espacial" #Calculando Indice de Moran para a variável resposta
moran(dados2$y, W.list, n = nrow(dados2), S0 = Szero(W.list))$I## [1] 0.07254203 #H0: nao ha dependencia espacial
#H1: ha dependencia espacial
IM=moran.test(dados2$y, W.list)
#IM$estimate[1]
IM##
## Moran I test under randomisation
##
## data: dados2$y
## weights: W.list
##
## Moran I statistic standard deviate = 1.2894, p-value = 0.09863
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.072542027 -0.010989011 0.004197009 alfa = 0.05
pvalor = IM$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("Não rejeito a hipotese de independencia espacial.")
}## [1] "Não rejeito a hipotese de independencia espacial." #######################################################################
# Ajustando o Modelo Espacial proposto por Leroux (CAR Leroux)
# sob o ponto de vista da Inferencia Bayesiana
# Y_i ~ N(X_i %*% beta + phi, tau2)
# c = (rho * sum_j=1^n(w_ij) + 1 - rho)
# phi_i|phi_i-1 ~ N(rho * sum_j=1^n(w_ij * phi_j) / c, nu2 / c)
#######################################################################
#Definindo as variaveis necessarias para o MCMC
n.aquec <- 100 #numero de iteracoes para o aquecimento
n.espac <- 150 #numero do espacamento
n.post <- 1000 #tamanho da amostra
n.ite <- n.aquec + n.post*n.espac #total de iteracoes
#Definindo os hiperparametros da distribuicao a priori
a_tau = 2
b_tau = 1
V_beta = 1000
a_V = 2
b_V = 1
# Ajustando o modelo via CARBayes
ajuste_CAR <- S.CARleroux(formula=dados2$y ~ dados2$x-1, family="gaussian",
W=W, burnin=n.aquec, n.sample=n.ite,
thin = n.espac, prior.mean.beta=rep(0,K),
prior.var.beta=rep(V_beta,K),
prior.tau2 = c(a_tau, b_tau),
prior.nu2 = c(a_V, b_V), verbose=FALSE)
cbind( ajuste_CAR$summary.results,
c(beta, V, tau2, rho) )## Mean 2.5% 97.5% n.sample % accept n.effective Geweke.diag
## dados2$x 1.9386 1.7416 2.1422 1000 100.0 1000.0 -0.1 2.0
## nu2 0.8092 0.5890 1.0807 1000 100.0 1000.0 0.4 1.0
## tau2 0.3596 0.1602 0.6672 1000 100.0 1000.0 -0.6 1.0
## rho 0.8998 0.6504 0.9944 1000 90.8 661.8 0.3 0.9 residuos_CAR = unname(ajuste_CAR$residuals$response)
#Erro médio quadrático (EQM)
mean((ajuste_CAR$fitted.values-y)^2)## [1] 0.6794013 #Deviance Information Criterion (DIC)
unname(ajuste_CAR$modelfit["DIC"])## [1] 257.1104 #------------------------------------
# Alocando as amostras a posteriori
rho_post = ajuste_CAR$samples$rho
tau2_post = ajuste_CAR$samples$tau2
V_post = ajuste_CAR$samples$nu2
beta_post = ajuste_CAR$samples$beta
phi_post = ajuste_CAR$samples$phi
# Analisando convergencia
# Traços a posteriori
par(mar=c(4,6,4,4),mfrow=c(2,2))
traceplot(beta_post[,1],ylab=expression(beta[1]),xlab="iterações",
cex.lab=1.8,cex.axis=1.6)
traceplot(V_post,xlab="iterações",ylab=expression(V),
cex.lab=1.8,cex.axis=1.6)
traceplot(tau2_post,xlab="iterações",ylab=expression(tau[2]),
cex.lab=1.8,cex.axis=1.6)
traceplot(rho_post,xlab="iterações",
ylab=expression(rho),cex.lab=1.8,cex.axis=1.6)
# Autocorrelação
par(mar=c(5,5,2,2),mfrow=c(2,2))
acf(beta_post[,1],bty="n",ylab=expression(beta[1]),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(V_post,bty="n",ylab=expression(V),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(tau2_post,bty="n",ylab=expression(tau^2),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(rho_post,bty="n",ylab=expression(rho),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
# Histograma
par(mfrow=c(2,2))
hist(beta_post[,1],freq=F,bty="n",ylab="densidade",xlab= expression(beta[1]),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=beta[1],lwd=2,col="red")
abline(v=mean(beta_post[,1]),lwd=2,col="blue")
abline(v=quantile(beta_post[,1],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(beta_post[,1],0.975),lwd=2,lty=3,col="blue")
hist(tau2_post,freq=F,bty="n",ylab="densidade",xlab=expression(tau^2),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=tau2,lwd=2,col="red")
abline(v=mean(tau2_post),lwd=2,col="blue")
abline(v=quantile(tau2_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(tau2_post,0.975),lwd=2,lty=3,col="blue")
hist(V_post,freq=F,bty="n",ylab="densidade",xlab="V",
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=V,lwd=2,col="red")
abline(v=mean(V_post),lwd=2,col="blue")
abline(v=quantile(V_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(V_post,0.975),lwd=2,lty=3,col="blue")
hist(rho_post,freq=F,bty="n",ylab="densidade",xlab=expression(phi),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=V,lwd=2,col="red")
abline(v=mean(rho_post),lwd=2,col="blue")
abline(v=quantile(rho_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(rho_post,0.975),lwd=2,lty=3,col="blue")
# Analisando os efeitos espaciais
phi.names=c(expression(phi[1]),expression(phi[2]),expression(phi[3]),
expression(phi[4]),expression(phi[5]),expression(phi[6]),
expression(phi[7]),expression(phi[8]),expression(phi[9]),
expression(phi[10]),expression(phi[11]),expression(phi[12]),
expression(phi[13]),expression(phi[14]),expression(phi[15]),
expression(phi[16]),expression(phi[17]),expression(phi[18]),
expression(phi[19]),expression(phi[20]),expression(phi[21]),
expression(phi[22]),expression(phi[23]),expression(phi[24]),
expression(phi[25]),expression(phi[26]),expression(phi[27])
)
# Traços de 16 phi
seqs = sort(sample(1:n, 16))
par(mar=c(4,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
traceplot(phi_post[,seqs[k]],xlab="iterações",ylab=phi.names[seqs[k]],
cex.lab=1.8,cex.axis=1.6)
abline(h=phi[seqs[k]],lty=1,lwd=2,col="red")
}
# Gráficos de autocorrelaçãoo dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(2,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
acf(phi_post[,seqs[k]],bty="n",ylab=phi.names[seqs[k]],main="",xlab="defasagem",cex.lab=1.8,cex.axis=1.6)
}
# Histograma da distribuição a posteriori dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(2,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
hist(phi_post[,k],freq=F,bty="n",ylab="densidade",xlab=phi.names[seqs[k]],main="",cex.lab=1.8,cex.axis=1.6)
abline(v=mean(phi_post[,k]),lwd=2,col="blue")
abline(v=quantile(phi_post[,k],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(phi_post[,k],0.975),lwd=2,lty=3,col="blue")
abline(v=phi[k],lwd=2,col="red")
}
# Estimativas Intervalares e Pontuais dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(4,6,4,4),mfrow=c(1,1))
media.phi = NULL
q.phi = matrix(NA,n,3)
for(i in 1:n){
media.phi[i] = mean(phi_post[,i])
q.phi[i,] = quantile(phi_post[,i],c(0.025,0.5,0.975))
}
plot(media.phi, ylim=c(min(q.phi),max(q.phi)),lwd=2,bty="n",axes=F,xlab="",ylab=expression(phi),cex.lab=1.8,cex.axis=1.6)
axis(2,cex.lab=1.4,cex.axis=1.4)
axis(1,at=1:n, label=1:n, cex.lab=1.4,cex.axis=1.4,las=2)
for(j in 1:n){
segments(x0=j,x1=j,y0=q.phi[j,1],y1=q.phi[j,3],lwd=2)
}
abline(h=0, lwd=2, col="red")
points(phi, pch="x", lwd=2)
# analisando os residuos do modelo espacial
# Histograma e qqplot dos res?duos do modelo espacial
par(mar=c(6,5,6,2),mfrow=c(1,2))
hist(ajuste_CAR$residuals$response,freq=F,cex.lab=1.8,cex.axis=1.6,ylab="densidade",xlab="resíduos",main="")
qqnorm(ajuste_CAR$residuals$response,cex.lab=1.4,cex.axis=1.4,ylab="quantis amostrais",xlab="quantis teóricos",main="",lwd=2,bty="n")
qqline(ajuste_CAR$residuals$response,col="red",lwd=2)
#analisando normalidade dos residuos
shapiro.test(ajuste_CAR$residuals$response)##
## Shapiro-Wilk normality test
##
## data: ajuste_CAR$residuals$response
## W = 0.97745, p-value = 0.1114 ks.test(ajuste_CAR$residuals$response, "pnorm",
mean(ajuste_CAR$residuals$response), sd(ajuste_CAR$residuals$response))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: ajuste_CAR$residuals$response
## D = 0.080119, p-value = 0.5685
## alternative hypothesis: two-sided #Calculando Indice de Moran para os resíduos
moran(ajuste_CAR$residuals$response, W.list, n = nrow(dados2), S0 = Szero(W.list))$I## [1] 0.2840103 #H0: nao ha dependencia espacial
#H1: ha dependencia espacial
IM=moran.test(ajuste_CAR$residuals$response, W.list)
IM##
## Moran I test under randomisation
##
## data: ajuste_CAR$residuals$response
## weights: W.list
##
## Moran I statistic standard deviate = 4.5697, p-value = 2.442e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.284010266 -0.010989011 0.004167428 alfa = 0.05
pvalor = IM$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("Não rejeito a hipotese de independencia espacial.")
}## [1] "Rejeito a hipotese de independencia espacial"2.3 Modelo de Regressão Poisson
Considere o seguinte modelo: \[\begin{eqnarray} Y_i | \mu_i &\sim& Pois(\mu_i), i=1, \ldots, n\\ g(\mu_i) &=& \log(\mu_i) = \boldsymbol{x}_i^T \boldsymbol{\beta} + O_i\\ \boldsymbol{\beta} &\sim& N_K(\boldsymbol{\mu}_\beta, \boldsymbol{\Sigma}_\beta) \end{eqnarray}\] sendo
\(\boldsymbol{x}_i\) o vetor de covariáveis da \(i\)-ésima região,
\(\boldsymbol{\beta}\) o vetor de efeitos das covariáveis e
\(O_i\) um parâmetro offset.
O termo offset pode ser entendido como uma reescalação nos dados.
Por exemplo, suponha que \(Y_i\) representa o número de casos de uma certa doença na região \(i\) e que as regiões tenham tamanhos populacionais diferentes. Encontrar 10 casos em uma região que tem 100 pessoas é diferente de encontrar 10 casos em uma região com 1.000 pessoas. Sendo assim, compensamos esse problema incluindo o termo \(O_i\) da seguinte forma no modelo apresentado acima: \[O_i = \log(\mbox{População total na região i}).\]
- Pode se incluir uma coleção de covariáveis espaciais neste modelo. Desta forma, o modelo não foi estruturado para ter dependência espacial mas dependendo das covariáveis pode capturar algum comportamento espacial.- Mais informações sobre este modelo podem ser obtidas em Agresti (2015).
#####################################################
#Simulando um valor observado para cada uma das regiões
#sob o modelo Poisson com distribuição a priori CAR Leroux
#####################################################
dados3 = Regiao
#Gerando um conjunto de dados
n = nrow(W)
beta = 2 #efeito da covariavel
K = length(beta)
x = rnorm(n, 0, sqrt(1))
I = diag(n) #matriz identidade
tau2 = 1 #variancia do termo espacial
Nviz = apply(W,1,sum) #total de vizinhos da regiao i
D = diag(Nviz)
rho = 0.9 #parametro de correlacao espacial
Sigma = tau2*solve(rho * (D-W) + (1-rho)*I) #matriz de covariancia do efeito espacial
phi = mvrnorm(1, rep(0,n), Sigma) #efeito espacial
mu = exp( x * beta + phi )
y = rpois(n, mu)
#plotando a variavel resposta
par(mfrow=c(1,1))
barplot(table(y), main="", cex.lab=1.4, cex.axis=1.4, ylab="frequência",
xlab="y")
dados3$y = y
dados3$x = x
#comparando a variável resposta e a covariável
#plotando a variável resposta
g1=ggplot() +
geom_sf(data = dados3, aes(fill=y), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
#plotando a covariável
g2=ggplot() +
geom_sf(data = dados3, aes(fill=x), color= NA, size=.15) +
labs(subtitle="", size=8) +
theme_minimal() +
no_axis
gridExtra::grid.arrange(g1, g2, top = '',ncol=2)
#Analisando phi
hist(phi, freq=F, main="Histograma", xlab=expression(phi),
ylab="densidade", cex.lab=1.4, cex.axis=1.4)
#atraves do Indice de Moran
alfa=0.05
pvalor = moran.test(phi, W.list)$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("Não rejeito a hipotese de independencia espacial.")
}## [1] "Rejeito a hipotese de independencia espacial"#Calculando Indice de Moran para a variável resposta
moran(dados3$y, W.list, n = nrow(dados3), S0 = Szero(W.list))$I## [1] -0.01484287#H0: nao ha dependencia espacial
#H1: ha dependencia espacial
IM=moran.test(dados3$y, W.list)
#IM$estimate[1]
IM##
## Moran I test under randomisation
##
## data: dados3$y
## weights: W.list
##
## Moran I statistic standard deviate = -0.099928, p-value = 0.5398
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.014842871 -0.010989011 0.001487373alfa = 0.05
pvalor = IM$p.value
if(pvalor<alfa){print("Rejeito a hipotese de independencia espacial")}else{
print("Não rejeito a hipotese de independencia espacial.")
}## [1] "Não rejeito a hipotese de independencia espacial."#######################################################################
# Ajustando o Modelo Espacial proposto por Leroux (CAR Leroux)
# sob o ponto de vista da Inferencia Bayesiana
# Y_i ~ Pois(mu_i)
#log(mu_i) = x_i beta + phi_i
# c = (rho * sum_j=1^n(w_ij) + 1 - rho)
# phi_i|phi_i-1 ~ N(rho * sum_j=1^n(w_ij * phi_j) / c, nu2 / c)
#######################################################################
#Definindo as variaveis necessarias para o MCMC
n.aquec <- 100 #numero de iteracoes para o aquecimento
n.espac <- 150 #numero do espacamento
n.post <- 1000 #tamanho da amostra
n.ite <- n.aquec + n.post*n.espac #total de iteracoes
#Definindo os hiperparametros da distribuicao a priori
a_tau = 2
b_tau = 1
V_beta = 1000
# Ajustando o modelo via CARBayes
ajuste_P_CAR = S.CARleroux(formula=dados3$y ~ dados3$x-1,
family="poisson",
W=W, burnin=n.aquec, n.sample=n.ite,
thin = n.espac, prior.mean.beta=rep(0,K),
prior.var.beta=rep(V_beta,K),
prior.tau2 = c(a_tau, b_tau), verbose=FALSE)
cbind( c(beta, tau2, rho), ajuste_P_CAR$summary.results)## Mean 2.5% 97.5% n.sample % accept n.effective Geweke.diag
## dados3$x 2.0 1.5532 1.3128 1.7665 1000 99.9 392.3 -0.1
## tau2 1.0 0.8691 0.4955 1.4484 1000 100.0 800.7 -2.9
## rho 0.9 0.6704 0.2717 0.9612 1000 94.7 319.5 0.9residuos_CAR = unname(ajuste_P_CAR$residuals$response)
#Erro médio quadrático (EQM)
mean((ajuste_P_CAR$fitted.values-y)^2)## [1] 0.9419514#Deviance Information Criterion (DIC)
unname(ajuste_P_CAR$modelfit["DIC"])## [1] 249.7501#------------------------------------
# Alocando as amostras a posteriori
rho_post = ajuste_P_CAR$samples$rho
tau2_post = ajuste_P_CAR$samples$tau2
V_post = ajuste_P_CAR$samples$nu2
beta_post = ajuste_P_CAR$samples$beta
phi_post = ajuste_P_CAR$samples$phi
# Analisando convergencia
# Traços a posteriori
par(mar=c(4,6,4,4),mfrow=c(2,3))
traceplot(beta_post[,1],ylab=expression(beta[1]),xlab="iterações",
cex.lab=1.8,cex.axis=1.6)
traceplot(tau2_post,xlab="iterações",ylab=expression(tau[2]),
cex.lab=1.8,cex.axis=1.6)
traceplot(rho_post,xlab="iterações",
ylab=expression(rho),cex.lab=1.8,cex.axis=1.6)
# Autocorrelação
acf(beta_post[,1],bty="n",ylab=expression(beta[1]),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(tau2_post,bty="n",ylab=expression(tau^2),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
acf(rho_post,bty="n",ylab=expression(rho),main="",xlab="defasagem",
cex.lab=1.8,cex.axis=1.6)
# Histograma
par(mfrow=c(1,3))
hist(beta_post[,1],freq=F,bty="n",ylab="densidade",xlab= expression(beta[1]),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=beta[1],lwd=2,col="red")
abline(v=mean(beta_post[,1]),lwd=2,col="blue")
abline(v=quantile(beta_post[,1],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(beta_post[,1],0.975),lwd=2,lty=3,col="blue")
hist(tau2_post,freq=F,bty="n",ylab="densidade",xlab=expression(tau^2),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=tau2,lwd=2,col="red")
abline(v=mean(tau2_post),lwd=2,col="blue")
abline(v=quantile(tau2_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(tau2_post,0.975),lwd=2,lty=3,col="blue")
hist(rho_post,freq=F,bty="n",ylab="densidade",xlab=expression(phi),
main="",cex.lab=1.8,cex.axis=1.6)
abline(v=V,lwd=2,col="red")
abline(v=mean(rho_post),lwd=2,col="blue")
abline(v=quantile(rho_post,0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(rho_post,0.975),lwd=2,lty=3,col="blue")
# Analisando os efeitos espaciais
phi.names=c(expression(phi[1]),expression(phi[2]),expression(phi[3]),
expression(phi[4]),expression(phi[5]),expression(phi[6]),
expression(phi[7]),expression(phi[8]),expression(phi[9]),
expression(phi[10]),expression(phi[11]),expression(phi[12]),
expression(phi[13]),expression(phi[14]),expression(phi[15]),
expression(phi[16]),expression(phi[17]),expression(phi[18]),
expression(phi[19]),expression(phi[20]),expression(phi[21]),
expression(phi[22]),expression(phi[23]),expression(phi[24]),
expression(phi[25]),expression(phi[26]),expression(phi[27])
)
# Traços de 16 phi
seqs = sort(sample(1:n, 16))
par(mar=c(4,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
traceplot(phi_post[,seqs[k]],xlab="iterações",ylab=phi.names[seqs[k]],
cex.lab=1.8,cex.axis=1.6)
abline(h=phi[seqs[k]],lty=1,lwd=2,col="red")
}
# Gráficos de autocorrelaçãoo dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(2,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
acf(phi_post[,seqs[k]],bty="n",ylab=phi.names[seqs[k]],main="",xlab="defasagem",cex.lab=1.8,cex.axis=1.6)
}
# Histograma da distribuição a posteriori dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(2,6,4,2),mfrow=c(4,4))
for(k in 1:length(seqs))
{
hist(phi_post[,k],freq=F,bty="n",ylab="densidade",xlab=phi.names[seqs[k]],main="",cex.lab=1.8,cex.axis=1.6)
abline(v=mean(phi_post[,k]),lwd=2,col="blue")
abline(v=quantile(phi_post[,k],0.025),lwd=2,lty=3,col="blue")
abline(v=quantile(phi_post[,k],0.975),lwd=2,lty=3,col="blue")
abline(v=phi[k],lwd=2,col="red")
}
# Estimativas Intervalares e Pontuais dos efeitos espaciais sorteados aleatoriamente
# do modelo espacial
par(mar=c(4,6,4,4),mfrow=c(1,1))
media.phi = NULL
q.phi = matrix(NA,n,3)
for(i in 1:n){
media.phi[i] = mean(phi_post[,i])
q.phi[i,] = quantile(phi_post[,i],c(0.025,0.5,0.975))
}
plot(media.phi, ylim=c(min(q.phi),max(q.phi)),lwd=2,bty="n",axes=F,xlab="",ylab=expression(phi),cex.lab=1.8,cex.axis=1.6)
axis(2,cex.lab=1.4,cex.axis=1.4)
axis(1,at=1:n, label=1:n, cex.lab=1.4,cex.axis=1.4,las=2)
for(j in 1:n){
segments(x0=j,x1=j,y0=q.phi[j,1],y1=q.phi[j,3],lwd=2)
}
abline(h=0, lwd=2, col="red")
points(phi, pch="x", lwd=2)
Referências
Agresti, Alan. 2015. Foundations of Linear and Generalized Linear
Models. 1º ed. Wiley.
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.
Leroux, Brian, Xingye Lei, e Norman Breslow. 2000. "Estimation of
disease rates in small areas: A new mixed model for spatial
diference". Em Statistical Models in Epidemiology, the
Environment, and Clinical Trials, editado por M. Elizabeth Halloran
e Donald Berry, 179–91. Spinger New York.