# 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" 

# Primeira versão do laboratório desenvolvido por Virginia Ragoni (2009), 
# Revisão elaborada por Flávia Feitosa e alunos da turma SER-301 (2011) 

#========================================================================
#                         PASSOS PRELIMINARES
#========================================================================

# 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


#========================================================================
#        IMPORTAÇÃO DOS DADOS A SEREM UTILIZADOS NAS ANÁLISES
#========================================================================

# São apresentadas 3 opções: 
# (1) com pacote 'rgdal'; 
# (2) Com pacote 'aRT'; 
# (3) Manualmente, exportando texto da tabela e matriz de vizinhança  
#     a partir de softwares como o TerraView ou o GeoDa. 


# ================== OPÇÃO 1: "rgdal" ===================================

install.packages('rgdal')
require('rgdal')

# Criando um objeto "Spatial Polygon DataFrame" a partir do shapefile
# objeto <-readOGR('.','nomedacamada')
# objeto <-readOGR('c:\\teste\\pastadoshp\\','nomedacamada')

# Exemplo no Linux: 
# objeto <- readOGR('/home/mauricio/Desktop/INPE/analise_espacial/aulas/columbus/','columbus')


tcolumbus<-readOGR('.','columbus')

# Criando um objeto "Spatial Polygon"
spols<-polygons(tcolumbus)

# Criando um objeto "Dataframe"
columbus<-data.frame(tcolumbus)

# Explorando um pouco estes objetos
class(spols)
plot (spols)
spplot(tcolumbus)
spplot(tcolumbus[1])

# ================== OPÇÃO 2: "aRT" ===================================

# Instruções para instalar o aRT (e-mail do Pedro)
# (1)Instalar o R 2.13.1 (http://cran-r.c3sl.ufpr.br/bin/windows/base/R-2.13.1-win.exe) em uma pasta que você tenha permissão (por exemplo, Meus Documentos, d:, e:) - não pode ser c:/arquivos de programas.
# Instalar o TerraView 4.0.0 ou 4.1.0 (http://www.dpi.inpe.br/terraview)
# Instalar o pacote sp: Executando o R, clique em Packages, Install Packages, Brazil (qualquer um), sp, OK.
# Instalar o pacote aRT: Na linha de comando do R, digite o seguinte comando e depois pressione Enter:
# install.packages("aRT", contrib="http://www.leg.ufpr.br/aRT")
# Instalar o MySQL: MySQL Community Server 5.5.16, Windows (x86, 32-bit), MSI Installer, disponível em http://www.mysql.com/downloads/
# Após instalar o MySQL, configure o servidor com as opções padrão. Um vídeo de como instalar e configurar está disponível em
# http://www.youtube.com/watch?v=oylVvq-3Tfk
# Note que ele trabalha com configurações avançadas. Toda vez que ele falar "you don't have to do this", você não precisa fazer para o curso. Note 
# tambem que se a configuração funcionar de primeira, você não precisa configurar novamente da mesma maneira que no vídeo.


# 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")

# thcolumbus <- createTheme(lcolumbus)

# Recuperar Tabela do Layer "lcolumbus" 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)

# ================== OPÇÃO 3: "Manual" ===================================

# A idéia é exportar tabela e matriz de vizinhança do terraview (pode ser no GeoDa, ou outro software)
# para posteriormente importá-la no R
# O aRT/ rgdal elimina estes passos de "exportação/importação"

# Criar um banco de dados no TerraView e importar os dados "columbus"
# 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 
# install.packages('spdep')
# library(spdep)
# polgal<-read.gal("columbus_neigh.GAL")

#========================================================================
#             EXPLORAÇÃO DE OBJETOS, CRIAÇÃO DE MATRIZ DE PESOS
#========================================================================

# Instalar todos os pacotes que serão utilizados no tutorial

install.packages('spdep')    
install.packages('lmtest')    
install.packages('lawstat')
install.packages('spgwr')


# ===== 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)
plot(polgal,coordinates(spols))
plot(spols,add=T)

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)
print(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)
 

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)