Dados de Área (continuação 2)

1 Analisando a correlação entre as localizações

Faça o seguinte exercício:

  • Obtenha o shapefile dos municípios de Minas Gerais usando o pacote geobr

  • Gere um dado simulado para estes municípios considerando um modelo SAR.

  • Replique essa geração 1000 vezes.

  • Calcule a correlação amostral e analise isso.

#carregando os pacotes
library(tidyverse)
library(rgdal)
library(maptools)
library(sp)
library(spdep)
library(spatialreg)
library(readxl)
library(geobr)
library(RColorBrewer)
library(MASS)
library(gridExtra)
library(ggiraph)
#install.packages("ggcorrplot")
library(ggcorrplot)

#Lendo o shapefile de dados
cod_estado    = 31
localizacoes  = read_municipality(code_muni = cod_estado,year=2018,simplified = T, showProgress = F)

# Construindo a matriz de vizinhança 
nb  = poly2nb(localizacoes,queen=TRUE)
n   = nrow(localizacoes)

# Ponderando a matriz de vizinhança (linhas somando 1)
nbw   = nb2listw(nb,style="W")
Wp    = as(nbw, "CsparseMatrix") %>% 
  as.matrix()

#analisando possiveis valores de rho
autovalores = eigen(Wp)
c(1/min(autovalores$values),1/max(autovalores$values))
## [1] -1.282833  1.000000
#Gerando o dado simulado sob o modelo SAR conforme solicitado
x             = matrix(1,n,2)   #matriz com intercepto + covariavel
x[,2]         = rnorm(n)        #covariável
beta          = c(5,1)          #efeito do intercepto + covariavel
V             = 0.5             #variancia
rho           = 0.9             #parametro de correlacao
nrep          = 1000            #total de replicações
e             = matrix(NA, n, nrep)
y_SAR         = matrix(NA, n, nrep)
minversa      = solve(diag(n)-rho*Wp)
mxbeta        = minversa %*% (x%*% beta)

for(r in 1:nrep){
  e[,r]              = rnorm(n, mean = 0, sd = sqrt(V))
  y_SAR[,r]          = mxbeta + minversa %*% e[,r]
}

##Calculando a correlação amostral
#cor_amostral2 <- cor(t(y_SAR))
##outra forma de calcular a correlação amostral
 cor_amostral = matrix(NA,n,n)
 for(i in 1:n){
   for(j in 1:n){
     cor_amostral[i,j] = cor(y_SAR[i,],y_SAR[j,])
   }
 }
#comparando as 2 formas de calcular
#summary(round(as.vector(cor_amostral-cor_amostral2),5))

##Analisando as correlações calculadas
#http://www.sthda.com/english/wiki/ggcorrplot-visualization-of-a-correlation-matrix-using-ggplot2
#x11()
##Plotando apenas os valores abaixo da diagonal principal
ggcorrplot(cor_amostral, method = "circle", type = "lower", pch=1)

##Deixando em branco as correlações não significativas
p.mat <- cor_pmat(t(y_SAR))
ggcorrplot(cor_amostral, p.mat = p.mat, hc.order = TRUE,  type = "lower", insig = "blank", method = "circle", pch=1)

#outras análises
#round(cor_amostral[1,],3)
#which(cor_amostral[1,]!=0)
#which(Wp[1,]!=0)
#cor_amostral[1,which(Wp[1,]!=0)]
#summary(cor_amostral[1,-c(1,which(Wp[1,]!=0))])
#cor_amostral[1,which(cor_amostral[1,]>0.5)]

par(mfrow=c(1,2))
hist(cor_amostral, freq=F, main="", ylab="densidade")
hist(cor_amostral[Wp!=0], freq=F, main="", ylab="densidade")

mean(cor_amostral[Wp!=0])
## [1] 0.7021207