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