# Aula 7 - Modelos de regressao linear multipla

getwd()
rm(list=ls())

# carregando os pacotes necess?rios
library(lattice)
library(sf)
library(spdep)

## 1 - leitura do banco de dados - ja arrumado!

### variaveis selecionadas apos eliminacao da colinearidade e ja no mapa dos munic do estado de SP
# "VAR_15","VAR_16","VAR_26","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70"

mapa <- st_read("./bd_aula_7/mapa_ca_mama_esp.shp")
plot(st_geometry(mapa))
print(mapa)
str(mapa)

## 2 - modelo com todas as covar

mod1a <- lm(b_TX_PAD ~ b_VAR_15+b_VAR_16+b_VAR_26+b_VAR_41+b_VAR_42+b_VAR_46+b_VAR_66+b_RZ_VAR_6+b_VAR_70, mapa) # mod de regress linear multipla
summary(mod1a) #sum?rio do modelo

# 3 - Selecao do modelo de menor AIC

## partir deste modelo, sem colinear entre as var indep, para chegar no mod final

# Seleção do modelo com menor AIC.
# AIC - Akikle Information Criterium = - 2 * log veross + 2 * número de paremetros do modelo

mod2a <- step(mod1a, trace = 1) # trace = mostrar resultados conforme roda o comando!

# Step:  AIC=5186.34
# TX_PAD ~ VAR_15 + VAR_42 + VAR_66

# var do mod: 
# VAR_15: % pop de 18 a 24 sem instr ou ate ciclo funda incom
# VAR_42: taxa de desemp
# VAR_66: renda per capta

summary(mod2a)

# Discussao - as varias formas de modelar!!!

# 4 - Analise dos residuos do modelo de regressao obito (mod2a)

# 4.1 Normalidade dos res

plot(mod2a, 2) # QQ-plot - avaliar normalidade dos residuos
ks.test(mod2a$residuals, mean(mod2a$residuals), sd(mod2a$residuals)) # teste para a normalidade

# 4.2 Autocorrelacao dos residuos

# construcao da matriz de viz - contiguidade
mapa.nb <- poly2nb(mapa) # matriz de viz com 1s e 0s
summary(mapa.nb)
str(mapa.nb)
print(mapa.nb)

# matriz de pesos
# style = "W" - soma dos pesos por linha = 1
# zero.policy=T - permite soma zero - algumas feicoes sem vizinhos

mapa.nbw <- nb2listw(mapa.nb, style="W", zero.policy = T) # matriz de pesos;
summary(mapa.nbw,zero.policy=T)
print(mapa.nbw,zero.policy=T)
str(mapa.nbw)

# Calc do I de Moran dos residuos do modelo

lm.morantest(mod2a, mapa.nbw, zero.policy = T)
# residuos sem dep espacial I = 0.079, p = 0.0003 # residuos com dep espacial - problema!!!

# 4.3 Residuos vs var indep - Homocedasticidade

# plot dos riduos vs valores preditos

plot(mod2a, 1) # outliers???

E0<-resid(mod2a) # colocando os res em um vetor E0

# plot dos residuos vs var indepent q entrarm no modelo

str(mapa)

# VAR_15
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_15 ,y = E0, xlab = "% pop 18-24 instr < ciclo fund", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
#Ok

# VAR_42
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_42 ,y = E0, xlab = "Taxa de desemprego", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# ok!

# VAR_66
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_66 ,y = E0, xlab = "Renda per capta", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# ok???

# plot dos residuos vs var indepent q nao entrarm no modelo
# se aparecer padr, a var, mesmo q nao tenha entrado no mod, tera q ser, de alguma forma, considerada
# VAR_16, VAR_26, VAR_41, VAR_46, RZ_VAR_68 e VAR_70

# VAR_16
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_16 ,y = E0, xlab = "VAR_16", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# ok!

# VAR_26
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_26 ,y = E0, xlab = "VAR_26", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# ok!

# VAR_41
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_41 ,y = E0, xlab = "Raz?o de renda", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# parece q temos prob. aqui!

# VAR_46
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_46, y = E0, xlab = "Taxa de analfabetismo em >= 15 anos", ylab = "Residuals", cex.lab = 1.5)
abline(h = 0, lty = 2)
#Ok

# RZ_VAR_68 - PIB per capta
par(mfrow = c(1,1))
plot(x = mapa$b_RZ_VAR_6 ,y = E0, xlab = "RZ_VAR_68", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# este parece q n?o ficou bom!

# VAR_70
par(mfrow = c(1,1))
plot(x = mapa$b_VAR_70 ,y = E0, xlab = "VAR_70", ylab = "Residuals",cex.lab = 1.5)
abline(h = 0, lty = 2)
# ok????