### AULA 8 REGRESSÃO ESPACIAL NO R

# remover objetos anteriores
#rm(list=ls())

# setar diretório

getwd() # verificação do diretório setado

# carregar arq de objetos (caso necessário!)

#load("aula_8.Rdata")

# carregando pacotes

library(spdep)
library(spatialreg)

#### Banco de dados

# vamos trabalhar com bd 'columbus' com informações sobre crimes
# os detalhes deste bd estão no arquivo 'columbus.html'

data(columbus) # para carregar o bd (está disponível junto com o 'spdep')
print(columbus)
#summary(columbus)

# após carregar 'columbus', cinco objetos ficam disponíveis

#ls()

# 'columbus' - data frame com 22 variáveis

# 'col.gal.nb' - objeto da classe 'nb' - lista de vetores, um para cada unidade espacial
# contendo a numeração sequencial dos vizinhos (critério - 'queen contiguity')
# OBS.: este obj já estava pronto no R, caso contrário, teríamos q criá-lo no GeoDa

str(col.gal.nb)

# 'coords' - centróides das coord - uso na constr de matriz de pesos com base na dist

str(coords)

# 'polys' e 'bbs' - o arq de contorno - uso para criar inform sobre contiguidade

print(polys)# informações sobre cada polígino q compõem a área de estudo
print(bbs)# dist aos 4 viz mais prox

### 2 RODANDO UMA REGRESSÃO

### 2.1 REGRESSÃO OLS

# Rodar a regressão OLS, para o bd 'Columbus' - função 'lm'
# var dependente = CRIME
# var independentes: INC (incone) e HOVAL (home value)
# aqui, vamos considerar que a análise exploratória e seleção de variáveis já foram realizadas!

columbus.lm <- lm(CRIME ~ INC + HOVAL, data = columbus)
columbus.lm # aparece um resumo dos resultados
summary(columbus.lm) # para ver os resultados completos

# interpretação dos resultados obtidos

### 2.2 PESOS ESPACIAIS

# função 'lm.morantest' - realiza teste estat I de Moran - autocorr espac dos res da regress
# opera em um objeto 'lm' - requer q um arq de pesos espac seja especificado
# essa função requer um objeto 'listw' - q tb contém valores para os pesos espaciais
# É necess converter um objeto 'nb' para um objeto 'listw' - usando função 'nb2listw'

col.listw <- nb2listw(col.gal.nb)

# verifique as difer entre 'col.gal.nb' e 'col.listw'

# verificando a estr de 'col.gal.nb'

#str(col.gal.nb)

# 'col.gal.nb' - lista de vetores - 1º vetor - elem 2 e 3
# depois da lista de obj - inform sobre os atributos dos obj, rótulos das áreas, etc

# verficando a estr de 'col.listw'

#str(col.listw)

# 'col.listw' - além do conteúdo de 'col.gal.nb' - inclui valores espec dos pesos espac

### 2.3 - TESTE I DE MORAN - AUTOCORREL ESPACIAL DOS RESÍDUOS

# função 'lm.morantest' 
# necessita o resultado da regress OLS (obj 'columbus.lm')
# necessita da informações de contiguidade fornecidas por um obj 'listw'

col.moran <- lm.morantest(columbus.lm,col.listw)
summary(col.moran)
print(col.moran) # teste significante - resíduos com dep espacial

# 'col.moran' - obj da classe 'htest' 
# nível de sign informado para a alternativa 'one-sided' ('greater')
# para a alternativa 'two-sided' - valor de p seria o dobro do valor informado
# espec da hip alterm - como abaixo

col.moran2 <- lm.morantest(columbus.lm,col.listw,alternative="two.sided")
print(col.moran2)

### IMPORTANTE

# esse teste leva em conta que são resíduos de uma regressão
# o teste I de Moran usual não leva isso em conta

str(columbus.lm)
col.val.res <- columbus.lm$residuals
col.val.res
summary(col.val.res) # resíduos com média zero

# teste para verific compor normal da var depend - Kolmogorov-Smirnov

ks.test(col.val.res,"pnorm",mean=mean(col.val.res), sd=sqrt(var(col.val.res))) 

# resíduos com distr normal

### 2.4 TESTE ESTAT MULTIPLICADOR DE LAGRANGE PARA AUTOCORR ESPACIAL

# teste I de Moran - não proporciona inform sobre modelos alternativos mais apropriados
# Teste Multiplicador de Lagrange - decisão entrre Mod do Erro Esp ou da Defasagem Esp

# função 'lm.LMtests' - arg: resultado da regress OLS e um objeto 'listw'
# além disso - especificar um vetor caracter (default é 'LMerror) - como abaixo

col.lagr <- lm.LMtests(columbus.lm,col.listw,test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
print(col.lagr)

# Comparando-se LMerr e LMlag - ambos são significantes
# Comparando-se as suas formas robustas: RMLerr com p = 0,85 e RLMlag com p = 0,07
# Sugestão de que o mod da defagem esp - alternativa mais provável

### 2.4 ESTIMAÇÃO DE MÁX VEROSSIM PARA MOD DE DEFASAGEM ESPACIAL

# função 'lagsarlm()' - estimação por máx veross (ML) de um mod de defasagem espacial
# arg: fórmula de regress, bd e objeto com pesos espaciais do tipo 'listw'

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

## comando 'larsarlm' está no pacote 'spatialreg' e não mais no spdep

# analisando-se os resultados

# parâmetro espacial autorregressivo (Rho)
# altamente significantes no teste assintótico t (p=0,0008)
# tb significante no teste da razão de veross (LR) (p=0,0037)
# Teste LR:
# hip nula: regr OLS; 
# hip altern: mod defasagem espacial ou erro

### 2.5 ESTIMAÇÃO DE MÁX VEROSSIM PARA MOD DE ERRO ESPACIAL

# função 'errorsarlm' - estim por máx veross (ML) de um mod de erro espacial
# arg: fórmula de regress, bd e objeto com pesos es´pac do tipo 'listw'

columbus.err <- errorsarlm(CRIME ~ INC + HOVAL,data=columbus,col.listw)
summary(columbus.err)
str(columbus.err)

# comparação do mod de defasagem com o de erro espacial - máx verossim
# log-verossim do mod de defasagem = -183,17 > log-veross do mod de erro = -184,16
# melhor modelo é o de defasagem espacial - como previsto pelo Teste Múltiplo de Lagrange.

#save.image("aula_8.Rdata")