# Laboratório utilizando algumas funcionalidades de regressão espacial
# contidas no pacote "spdep" em R, devenvolvido por Roger Bivand

# O laboratório foi baseado no artigo de Luc Anselin
# "An Introduction to Spatial Regression Analysis in R" 

# Laboratório desenvolvido por Virginia Ragoni (2009)
# e revisado por Flávia Feitosa (2011)

#========================================================================
#                PREPARAÇÃO E EXPLORAÇÃO DE OBJETOS
#========================================================================

# Criar uma pasta de trabalho e defini-la como área de trabalho no R
# Download arquivo columbus.zip (Columbus neighborhood crime) 
# em http://geodacenter.asu.edu/sdata
# Descompactar e salvar arquivos na área de trabalho

# ==== GERAR BANCO DE DADOS NO TERRAVIEW (USANDO aRT)=============

require(aRT)
aRTsilent(FALSE)

# Lembre-se de colocar sua senha aqui para conexão com o BD:
con = openConn(u="root",pass="fff0111") 

if(any(showDbs(con)=="columbus")) deleteDb(con, "columbus", force=TRUE)
db=createDb(con, "columbus")

lcolumbus = createLayer(db, "columbus")

addShape(lcolumbus, tab="columbus", file="columbus.shp", ID="POLYID")

tcolumbus = createTheme(lcolumbus)

# Recuperar Tabela do Layer "columbus" e criar data frame
tcolumbus <- openTable(lcolumbus, "columbus")
columbus <- aRT::getData(tcolumbus)

# Criar um objeto "Spatial Polygon" a partir do Layer columbus
spols <- getPolygons(lcolumbus)

class(spols)
plot (spols)


# ====================== Opção "Sem aRT"  ========================
# A idéia é exportar tabela e matriz de vizinhança do terraview para
# posteriormente importá-la no R
# O aRT elimina estes passos de "exportação/importação"

# Abrir o banco de dados "columbus" no TerraView e visualizar o dado

# Exportar tabela columbus.txt para área de trabalho do R: 
# Clicar na tabela com botão direito > Tabela > Exportar Tabela como...
# > Salvar tabela como "columbus.txt", separada por vírgula,
# no diretório de trabalho do R

# Criar matriz de vizinhança (contiguidade) no TerraView 
# Análise > Matrix de Proximidade > Criar Matriz de Proximidade
# Salvar no diretório de trabalho do R

# No R: Importar tabela e criar um dataframe

# columbus <- read.csv("columbus.txt", header=TRUE)

# Converter a matriz de vizinhança (externa) em um objeto contendo uma lista de vizinhos
# Necessário instalar e carregar pacote spdep 
# library(spdep)
# polgal<-read.gal("columbus_neigh.GAL")



# ===== R: EXPLORANDO OS DADOS ======================================

# Visualizar o objeto "columbus"
columbus
edit(columbus)

# Confirmar se o objeto "columbus" é um dataframe
class(columbus)
is.data.frame(columbus)  

# Vamos explorar um pouco o banco:
# (1)Verifique o nr. de linhas/colunas/atributos:
nrow (columbus)
ncol (columbus)
length (columbus)

# (2) Nomes dos atributos:
names(columbus)

# (3) Olhe os atributos das 5 primeiras linhas
columbus[1:5,]
print(columbus[1:5,])

# (4) Estatística descritiva para cada variável, e somente para a 9a. variável (CRIME)
summary(columbus)
summary(columbus[9])

# (5) Fazer o histograma para a variável crime
hist(columbus[,9])


# ===== CRIAR LISTA DE VIZINHOS E DE PESOS ============

# Criar lista de vizinhos a partir do objeto "Spatial Polygon"
# Necessário instalar e carregar pacote spdep 

library(spdep)
polgal<-poly2nb(spols)
# polgal<-poly2nb(spols, queen = FALSE)

class(polgal)
# No pacote spdep, as relações de vizinhança entre n observações 
# são representadas por um objeto da classe "nb" - Neighbour list

# Explorar o objeto 
polgal
summary(polgal)
is.list(polgal)
length(polgal)
polgal[[1]]
polgal[[1]][2]
polgal[[49]][1]

# Criar uma lista de pesos a partir da lista de vizinhos:
# Comando nb2listw
W_polgal<-nb2listw(polgal)

# Verificar a lista de pesos criada
summary(W_polgal)
class(W_polgal)

# Verificar os pesos contidos no objeto
W_polgal$weights
W_polgal$weights[1]


#========================================================================
#                         REGRESSÃO LINEAR
#========================================================================

# Variável resposta é CRIME
# Variável explicativa é INC
# CRIME = roubos residenciais e de veículoS
# INC  = rendimento familiar
# HOVAL = valor dos imóveis


# plotando a relação entre as variáveis
plot(CRIME~INC, data=columbus)
columbus.lm.inc <-lm(CRIME~INC,data=columbus)
abline(columbus.lm.inc)
summary(columbus.lm.inc)

plot(CRIME~HOVAL, data=columbus)
abline(columbus.lm.hoval<-lm(CRIME~HOVAL,data=columbus))
summary(columbus.lm.hoval)

# chamar a função apresenta alguns resultados
# lm(CRIME ~ INC, data=columbus)

lm(CRIME ~ INC + HOVAL, data=columbus)


# guarda a função como objeto da classe
class(lm)
#columbus.lm <- lm(CRIME ~ INC , data=columbus)

columbus.lm <- lm(CRIME ~ INC + HOVAL , data=columbus)

class(columbus.lm)
is.list(columbus.lm)
summary(columbus.lm)
names(columbus.lm)

# verificar o significado de cada variável.
 
summary(columbus.lm)    # para resultados detalhados

# Plotar o objeto spols (spatial polygons - columbus)
plot(spols)

# Plotar os residuos do modelo usando 4 intervalos
brks <- round(fivenum(columbus.lm$res), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(columbus.lm$res,brks)])
title(main="Mapa dos resíduos da regressão linear multipla CRIME~INC + HOVAL")
legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)

# gráfico dos resíduos x cada variável explicativa
hist(columbus.lm$residuals)
hist(columbus$CRIME)

# Diagnósticos da regressão
plot(columbus.lm$residuals, columbus$INC)
title(main= "Mapa dos resíduos x INC")
plot(columbus.lm$residuals, columbus$HOVAL)
title(main= "Mapa dos resíduos x HOVAL")


plot(columbus.lm$fit, columbus.lm$residuals)
title(main= "Gráfico dos valores ajustados X resíduos")
qqnorm(columbus.lm$residuals, ylab= "Residuos")
qqline(columbus.lm$residuals)


# alguns testes:

# de normalidade dos dados - Bera_Jarque test para não normalidade (hipótese nula é que é normal)
# Instalar biblioteca lawstat
library(lawstat)

rjb.test(columbus.lm$residuals, option="JB", crit.values="chisq.approximation",N=0)

# teste de Breusch-Pagan para heterocedasticidade (H0=homocedasticidade)
# carregar biblioteca lmtest
library (lmtest)
bptest(CRIME ~ INC + HOVAL, data=columbus)

# Em caso de heterocedasticidade, medidas de remediação podem incluir: 
# Tranformações das variáveis, ou
# Weighted Least Square Regression: http://itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm


#========================================================================
#                    TESTE DE MORAN PARA OS RESIDUOS
#========================================================================

#   O teste de Moran para autocorrelação espacial é aplicado aos resíduos da regressão.
#   a função é lm.morantest
#   um dos parâmetros dessa função é o arquivo de pesos espaciais.
#   esse arquivo é um objeto  listw.
#   listw não é o mesmo que o arquivo nb (col.gal.nb) - 
#   nb contém apenas os identificadores dos vizinhos e listw contem os pesos espaciais
#   é necessário transformar nb em listw (conforme feito em etapa anterior --> W_polgal)

tmoran <- lm.morantest (columbus.lm, W_polgal)
tmoran

# verifique o p-valor, tenha em mente de que este é um teste unilateral
# Para teste bilateral, faça:  

tmoran2 <- lm.morantest (columbus.lm, W_polgal, alternative="two.sided")

# Conclui-se que a hipótese nula de não correlação é rejeitada
# existe correlação espacial

##################################
# E agora?
#
# Utlizar os testes dos multiplicadores de lagrange que permitem distinguir os
# modelos spatial error e os modelos spatial lag
# 
# Ambos os testes e suas fomas robustas estão incluídas na função lm> LMtests.
# os testes devem ser definidos, caso contrário apenas o LMerror é calculado (default)
#################################

columbus.lagrange <- lm.LMtests(columbus.lm,W_polgal,test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
columbus.lagrange

# Primeiro deve-se observar qual dos LMerr e LMlag são significantes (p-valor).
# os dois são significantes. Quando isso ocorre deve-se comparar a versão robusta
# nesse caso, RMerr tem p-valor não significante (0.834) e RMlag é significante
# (0.05326) 
# O modelo a ser escolhido é o spatial lag 

#========================================================================
#       MODELO ESPACIAL DO TIPO SAR (Spatial Autoregressive Model)
#                        Spatial Lag Model
#========================================================================

# A estimação do modelo será por Maxima Verosemelhança (Maximum likelihood) 
#   y = rho W y + X beta + e


columbus.lag <- lagsarlm(CRIME~INC+HOVAL,data=columbus,listw=W_polgal)
summary(columbus.lag)
names(columbus.lag)

# o parâmetro RHO é significante - 
# o p-valor do teste-t assintótico (baseado na variância assintótica) é 0.0004
# O teste da razão de máxima verossimilhança (LR) é também significante 0.0022

# plotando o mapa de resíduos
plot(spols)
brks <- round(fivenum(columbus.lag$residuals), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(columbus.lag$residuals,brks)])
title(main="Mapa dos resíduos do modelo SAR : CRIME~INC+HOVAL")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)


#========================================================================
#       MODELO ESPACIAL DO TIPO CAR (Conditional Autoregressive Model)
#                        Spatial Error Model
#========================================================================

# A estimação do modelo sera por Maxima Verosemelhança (Maximum likelihood) 
# y = X ß + u
# u = Lambda W u + epsilon

columbus.err<-errorsarlm(CRIME~INC+HOVAL,data=columbus,listw=W_polgal)  
summary(columbus.err)
names(columbus.err)
class(columbus.err)

# o parâmetro LAMBDA é significante - 
# o p-valor do teste-t assintótico (baseado na variância assintótica) é 0.0002
# O teste da razão de máxima verossimilhança (LR) é também significante 0.01 

# Uma comparação entre os dois modelos deve ser baseada na 
# maximização do log da verossimilhança  e não no valor de R2 .
# O valor do log da verossimilhança para o modelo CAR é -184.16 e para o SAR
# é -182.67. O modelo CAR então é pior do que o SAR.

# O valor de AIC para o modelo CAR é 378.31 Para o SAR é 375.35.
# Esses valores não devem ser comparados entre si mas sim com a regressão simples
# cujo AIC é pior 382.75

 # plotando o mapa de resíduos

brks <- round(fivenum(columbus.err$residuals), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(columbus.err$residuals,brks)])
title(main="Mapa dos resíduos do modelo CAR : CRIME~INC+HOVAL")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)



#========================================================================
#          MODELO GWR (Geographically Weighted Regression)
#========================================================================

# Instalar e carregar pacote spgwr
 library(spgwr)

# seguindo o livro do Bivand

# calculando a largura banda fixa - para filtro adaptativo adapt=TRUE)
bw <- gwr.sel(CRIME~INC+HOVAL,data=columbus,coords=cbind(columbus$X, columbus$Y))
#bwa <- gwr.sel(CRIME~INC+HOVAL,data=columbus,coords=cbind(columbus$X, columbus$Y), adapt=TRUE)


gwr_columbus<-gwr(CRIME~INC+HOVAL,data=columbus,coords=cbind(columbus$X, columbus$Y), bandwidth= bw, gweight=gwr.Gauss, hatmatrix=TRUE)
#gwra_columbus<-gwr(CRIME~INC+HOVAL,data=columbus,coords=cbind(columbus$X, columbus$Y), bandwidth= bwa, gweight=gwr.Gauss, hatmatrix=TRUE)
 

#pairs(as(gwr_columbus$SDF, "data.frame")[,2:8], pch=".") # difícil de ver

  col_intercept= gwr_columbus$SDF[[2]]
  col_income = gwr_columbus$SDF[[3]]
  col_housing =  gwr_columbus$SDF[[4]]
  
hist(col_intercept)
hist(col_income)
hist(col_housing)

# mapeando os coeficientes da regressão para cada ponto
# intercepto
brks <- round(fivenum(col_intercept), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(col_intercept,brks)])
title(main="Mapa do intercepto")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)

#income
brks <- round(fivenum(col_income), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(col_income,brks)])
title(main="Mapa do coeficiente de INCOME")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)

#housing
brks <- round(fivenum(col_housing), digits=2)
brks
cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(col_housing,brks)])
title(main="Mapa do coeficiente de HOUSING")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)


 #resíduos?
 
col_res =  gwr_columbus$SDF[[8]]

 brks <- round(fivenum(col_res), digits=2)
 brks
 #cols <- grey(5:2/6)

 cols <- rev(heat.colors(4))
plot(spols,col = cols[findInterval(col_res,brks)])
title(main="Mapa dos resíduos")

legend(c(6, 8), c(13,15), legend = leglabs(brks,"<", ">="), fill = cols, bty = "n", cex = 0.9, y.intersp = 0.9)