# Aula 6 - Analise Exploratoria de Dados

# Aula elaborado com base no artigo de Zuur AF, Ieno EN, Elphick CS. A protocol for data exploration to avoid common statistical problems. Methods in Ecology & Evolution 2010, 1, 3-14.

getwd()
rm(list=ls())

#load(aula6.Rdata)

# instalar pacotes necessarios

# carregar pacotes

library(lattice)
#library(rgdal)
library(spdep)
library(car) # VIF
library(ggplot2)
library(caret)
library(gstat)
library(sf)
library(nlme)

# funcao desenvolvida por Zuur et al para ajudar na an explor 
# ficar junto do script
source("HighstatLibV9.R") 

# leitura do banco de dados - obitos por ca de mama munic do ESP
# artigo Diniz et al 2017

banco <- read.csv('./bd_aula6/banco.csv', sep=";",dec = '.')
head(banco)

# 1. Outliers em X e em Y

# Outlier: observ que tem um valor relativamente grande ou pequeno comparado com a maioria das observacoes
# outliers podem causar superdispersao (overfitting) em Modelos Lineares Generalizados
# Ferramentas utilizadas: 

# Boxplot - nem sempre os valores acima (ou abaixo) das linhas sao outliers

# Cleveland dotplot: 
# o numero da linha de uma observacao é plotado vs. o valor da observacao
# proporciona inform muito mais detalhadas do q o boxplot
# pontos q ficam no lado direito ou esquerdo do grafico podem corresp a outliers
# sao valores que sao considerados maiores ou menores do que o maioria das observacoes

# sao pontos que requerem investigacao
# checagem dos dados brutos
# verificacao sobre a coerencia e a razoabilidade do dado


## 1.1 Outliers na vari?vel resposta.

boxplot(banco$TX_PAD)
dotchart(banco$TX_PAD)

# comparando os dois graficos
#valores apontados como outliers pelo boxplot nao sao tao extremas qto pareceriam qdo olhamos o dotplot

## 1.2 Outliers nas variaveis independentes
str(banco)

MyVar <- c("VAR_15","VAR_16","VAR_17","VAR_25","VAR_26","VAR_27","VAR_31","VAR_41")
Mydotplot(banco[,MyVar]) # usando a funcao 'HighstatLibV9.R'
boxplot(banco$VAR_31) # valores acima da linha super do boxplot, neste caso, nao sao outliers


MyVar <- c("VAR_42","VAR_46","VAR_66","VAR_68","VAR_69","VAR_70")
Mydotplot(banco[,MyVar])

# temos outlier na VAR_42 e na VAR_68 - vamos verificar

banco[order(banco$VAR_42,decreasing=T),] # ordenacao
# VAR_42 - Taxa de desemprego = 36.52 - erro de digitacao!
# Qdo explicao mais provavel eh erro de medida - valores devem ser excluidos
# A sua presenca pode dominar a analise
## usar o View ou clicar no objeto em "Environment"

# corrigindo valor- correto é 6.52%

banco_nv <- banco
banco_nv['47','VAR_42']
banco_nv['47','VAR_42'] <- 6.52
banco_nv['47','VAR_42']

MyVar <- c("VAR_42","VAR_46","VAR_66","VAR_68","VAR_69","VAR_70")
Mydotplot(banco_nv[,MyVar])

banco[order(banco$VAR_68,decreasing=T),] # ordenacao
# VAR_68 - PIB per capta - eh o valor de Louveira - tem um PIB mais alto mesmo! Nao eh erro de digit.
# Opcao 1: fazer analise sem ele, mas considera-lo na apresentacao dos resultados (val pred, por ex.)
# Opcao 2: transformar a covariavel usando, por exemplo, raiz quadrada (o q foi feito no artigo)

# criando a covar RZ_VAR_68

banco_nv$RZ_VAR_68 <- sqrt(banco$VAR_68)
head(banco_nv)

MyVar <- c("VAR_42","RZ_VAR_68")
Mydotplot(banco_nv[,MyVar])

# exportando o banco_nv, caso seja necessario

write.csv(banco_nv,file="./bd_aula5/banco_nv.csv")

# 2. Homogeneidade da variancia

# importante em Analise de Variancia - ANOVA - comparacao de medias

# Analise de regressao: 
# comportamento dos residuos versus valores ajustados
# comportamento dos residuos versus valores da var indep - sem padrao
# possiveis solucoes: transf da var resp, minmos quadr generalizados

# 3. Normalidade

# algumas tecnicas estatisticas - normalidade da var dependente
# Regressao linear - normalidade dos residuos

# ferramentas: histogramas e QQ-plots

hist(banco_nv$TX_PAD)
qqnorm(banco_nv$TX_PAD)

# 4. Excesso de zeros

# Altas proporcoes de zero na var dependente
# producao de estimativas de parametros e de erros padroes enviesados
# considerar os Modelos Inflados de Zero (ZIM - Zero Inflated Models)

sum(banco_nv$TX_PAD == 0) / nrow(banco_nv) * 100 # 10%?de valores iguais a zero

# 5. Colinearidade entre as variaveis independentes

# refere-se ao prob de existencia de correlacao entre as var. indep.
# Se o prob for ignorado - produz modelos ruins
# exemplo - Tabela 1 do art de Zuur et al - pag 9

# Ferramentas

str(banco_nv)
pairs(banco_nv[,c(5:15,17,18,20)])
cor <- cor(banco_nv[,c(5:15,17,18,20)]) # criterio de corte - coef de correl >= 0.5 a 0.8

# exportanto matriz de correlacao

write.table(cor, file = "cor.txt", sep = ",", col.names = NA, qmethod = "double")

# outro modo usando o comando 'MYpairs" da funcao 'HighstatLibV9.R'

MyVar <- c("VAR_15","VAR_16","VAR_17","VAR_25","VAR_26","VAR_27","VAR_31","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_69","VAR_70")
Mypairs(banco_nv[,MyVar])

# VIF - variance inflation factor (pag 9, primeira coluna do artigo)
# criterio VIF < 3

# Remoção da colinearidade - utilizando o VIF

corvif(banco_nv[,MyVar])

# sai VAR_27

MyVar <- c("VAR_15","VAR_16","VAR_17","VAR_25","VAR_26","VAR_31","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_69","VAR_70")
corvif(banco_nv[,MyVar])

# sai VAR_69

MyVar <- c("VAR_15","VAR_16","VAR_17","VAR_25","VAR_26","VAR_31","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70")
corvif(banco_nv[,MyVar])

# sai VAR_17

MyVar <- c("VAR_15","VAR_16","VAR_25","VAR_26","VAR_31","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70")
corvif(banco_nv[,MyVar])

# sai VAR_25

MyVar <- c("VAR_15","VAR_16","VAR_26","VAR_31","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70")
corvif(banco_nv[,MyVar])

# entre var_31, % de pessoas com renda menor do q 1/2 sm, e var_66, renda per capta, ficamos com a segunda (mais representativa)

## podemos olhar a relação de cada uma delas com a var resposta, antes de decidir (ver a seguir)

# sai VAR_31

MyVar <- c("VAR_15","VAR_16","VAR_26","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70")
corvif(banco_nv[,MyVar]) ### OK! todas covar com VIF < 3


# 6. Relacao em var dependente (Y) e var indep (X)

# Plotagem da var dependente versus cd uma da var indep
# A ausencia de padroes claros entre X e cd uma da var indep nao significa q nao haja relacao entre elas.
# Signfica apenas q nao um ha claro relacionamento bivariado entre X e Y
# Modelagem com multiplas var explan pode proprocionar ainda um bom ajuste

### pode ajudar na decisao na escolha entre covariaveis colineares - entre duas ou mais, escolher a que tem "melhor" relacao com  variael dependente

# plotando Y em relacao a cd um dos Xs - usando o script 'HighstatLib.V9.R'

MyX  <- c("VAR_15","VAR_16","VAR_26","VAR_41", "VAR_42","VAR_46","VAR_66","RZ_VAR_68","VAR_70")
Myxyplot(banco_nv, MyX, "TX_PAD", MyYlab = "Tx Padron")

# VAR_46(analf) e VAR_66 (rend per cap) parecem ter relac com TX_PAD

# Outras utilidades import dos Scatter-plots:
 
# Visualizacao de relacionamento entre variaveis
# Detecta??o de observa??es q n?o est?o em conformidade com o padr?o geral entre duas vari?veis - outliers

# Qq observ q fique fora da n?vem preta de ptos necessita ser investigada
# podem ser dados de outra doen?a, erros de medida, erros de digita??o e tb podem ser valores corretos (como o PIB per cap de Louveira)
 
# 7. Existem interacoes?

# graficos do tipo 'coplot' sao excelentes ferramentas para visualizar a presenca potencial de interacoes

# vamos ler o arquivo dos retardos de diag e verificar se existe potencial interacao entre a diferenca da dist percorrida (Real - Ideal) 
# e a classificacao da TB (extra e pulm)

# isto eh, vamos verificar se a relacao entre o retardo do diag e a dif da dist percorrida muda em funcao da classif da TB

ret.diag <- read.table("./bd_aula6/ret_diag.txt", header=T)
head(ret.diag)
str(ret.diag)

coplot(RET_DIAG ~ DMAIOR_D_1 | factor(CLAS_P_EX), data = ret.diag,
       panel = function(x, y, ...) {
         tmp <- lm(y ~ x, na.action = na.omit)
         abline(tmp)
         points(x, y) })

# vemos q ha uma interacao em poten entre as duas var - as linhas nao sao paralelas

# para confirmar isso, devemos rodar um mod de regress?o e testar a exist de intera??o.

coplot(RET_DIAG ~ VZ_PR_SERV | factor(CLAS_P_EX), data = ret.diag,
       panel = function(x, y, ...) {
         tmp <- lm(y ~ x, na.action = na.omit)
         abline(tmp)
         points(x, y) })

# 8. As var categoricas estao balanceadas?

table(ret.diag$CLAS_P_EX) # mais ou menos
table(ret.diag$RELIG) # religiao 4 e 6 com uma unica obs cd - problema!!!
table(ret.diag$TIPO_ATENC)
table(ret.diag$ESCOL)

# 9. As obs da var depend sao independentes?

# pressuposto basico de muitas tec estat - obs da var Y devem ser indep uma das outras

# No banco de dados de casos de dengue de Varzea Paulista de 2007
# Avaliar se var dep tem dep espacial

# lendo o shape 'munic_esp_ob_pop_2010.shp'

mama <- st_read("./bd_aula6/munic_esp_ob_pop_2010.shp") # lendo mapa

class(mama)
plot(mama)
print(mama)
str(mama)
names(mama)

# construindo matriz de vizinhanca

mama.gal.q <- poly2nb(as(mama, "Spatial"),queen=TRUE)
mama.gal.q
summary(mama.gal.q)
card(mama.gal.q) # n? de viz de cd polig

mama.listw <- nb2listw(mama.gal.q, zero.policy=TRUE) # transf em pesos
summary(mama.listw)

?nb2listw

# moran da var depend

moran.test(mama$TX_PAD, mama.listw)

# na verdade, o q ira nos interessar eh se os resduos das modelagens tem ou nao depend espacial.

####################################################

save.image("aula6.Rdata")

#################################################