## Codigos apresentados na aula 4: Análise exploratória de dados ## ##################################################### ### Uso da função is.na e substituição de valores ### ##################################################### setwd("SEU DIRETÓRIO DE TRABALHO") getwd() ## Lendo a planilha com read.table aves.c <- read.table("aves_cerrado.csv", row.names=1, header=T, sep=";", dec=",", as.is=T) ## Ou com read.csv2 aves.c <- read.csv2("aves_cerrado.csv", row.names=1, as.is=T) ## Verificacao head(aves.c) tail(aves.c) str(aves.c) summary(aves.c) ## Agora vamos verificar os valores da coluna que será um fator table(aves.c$fisionomia) ## Há um erro de digitacao no codigo de fisionomia de uma parcela. ## Corrigindo: aves.c$fisionomia[aves.c$fisionomia=="ce"] <- "Ce" table(aves.c$fisionomia) class(aves.c$fisionomia) ## Convertendo para fator, que ordenamos da fisionomia mais aberta ## para a menos aves.c$fisionomia <- factor(aves.c$fisionomia, levels=c("CL","CC","Ce")) ## Verificando tudo novamente str(aves.c) summary(aves.c) ## VOLTAR P/ APRESENTACAO - primeiro problema: os NAs são de fato observações faltantes? ## Vamos ver para os urubus aves.c[aves.c$urubu==NA,] ## Ops, nao funciona ## O teste logico para NA é feito pela funcao is.na: is.na(aves.c) is.na(aves.c$urubu) ## Contando quantos NAs tem em cada coluna avesNAs <- is.na(aves.c[,2:4]) colSums(avesNAs) ## Quero saber qual é o registro de NA para urubu aves.c[is.na(aves.c$urubu),] ## Verificando a planilha descobrimos que na verdade são pontos ## em que não houve avistamento. Portanto, o valor correto é zero ## Quais sao os registros com este problema? ## E para NA em uma das tres aves aves.c[is.na(aves.c$urubu) | is.na(aves.c$carcara)| is.na(aves.c$seriema), ] ## Vamos guardar este pedaço da planilha num objeto temporario temp1 <- aves.c[is.na(aves.c$urubu)|is.na(aves.c$carcara)| is.na(aves.c$seriema),] ## E agora vamos corrigir estes valores aves.c$urubu[is.na(aves.c$urubu)] <- 0 ## Que é o mesmo que aves.c[is.na(aves.c$urubu), 2] <- 0 ## Ou ainda aves.c[is.na(aves.c[,2]), 2] <- 0 ## Continuando, para as outras aves aves.c$carcara[is.na(aves.c$carcara)] <- 0 aves.c$seriema[is.na(aves.c$seriema)] <- 0 ## Verificando aves.c[is.na(aves.c$urubu)|is.na(aves.c$carcara)| is.na(aves.c$seriema), ] ## Que devem ser agora valores com zero: aves.c[aves.c$urubu==0|aves.c$carcara==0|aves.c$seriema==0,] ## Comparando com o arquivo temporario criado acima temp1 # ok! ## Se eu quisesse saber o número de zeros aveszeros <- aves.c[,2:4]==0 aveszeros colSums(aveszeros) ## Uso da funcao na.omit ## Boa para funcoes que nao tem argumento na.rm (coisa <- c(1:5,NA,NA)) na.omit(coisa) mean(coisa) mean(coisa, na.rm=TRUE) mean(na.omit(coisa)) # se a funcao mean nao tivesse o argumento na.rm ########################################################### ### Resumo estatistico: medias, media truncada e mediana ### ############################################################ summary(aves.c[,2:4]) sd(aves.c[,2:4]) ## Funciona, mas está "condenada" (deprecated). A maneira recomendada é: sapply(aves.c[,2:4],sd) # aplica uma funcao a cada elemento de uma lista ## ou apply: apply(aves.c[,2:4], 2, sd) # data frames são listas, mas aceitam o # apply, por terem arranjo matricial ## Prosseguindo: mediana, media, media truncada apply(aves.c[,2:4], 2, median) apply(aves.c[,2:4], 2, mean) #ou a mesma coisa com o colMeans colMeans(aves.c[,2:4]) apply(aves.c[,2:4], 2, mean, trim=0.1) ## Quantis quantile(aves.c$urubu) ## O mesmo que o retornado pelo summary summary(aves.c$urubu) ## Mas vc pode mudar os quantis quantile(aves.c$urubu, probs= seq(from=0,to=1,by=0.1)) ############## OPCIONAL ######################## ### Parenteses: apply e sapply ### ## apply aplica a funcao FUN às margens de um objeto de arranjo matricial ## (classes matrix ou data frame) apply(aves.c[,2:4],1,range) # minimo e maximo por linhas apply(aves.c[,2:4],2,range) # minimo e maximo por colunas ## sapply: aplica uma funcao a cada elemento de uma lista ## Exeplo com uma lista com 4 resultados de regressoes ## Construindo a lista com os dados do objeto Anscombe asc.regr <- list(m1=lm(y1~x1, data=anscombe), m2=lm(y2~x2, data=anscombe), m3=lm(y3~x3, data=anscombe), m4=lm(y4~x4, data=anscombe)) names(asc.regr) names(asc.regr$m1) ## Coeficientes de cada regressao: intercepto e inclinacao ## (aqui com nome x1) sapply(asc.regr,coef) coef(asc.regr$m1) ############################ ### Analises univariadas ### ############################ ## Contagens de fatores: tabelas caixeta <- read.csv("caixeta.csv", as.is=T) names(caixeta) str(caixeta) table(caixeta$especie) sort(table(caixeta$especie), decreasing=T) table(caixeta$local) ## Graficos de barra para representar uma tabela barplot(sort(table(caixeta$especie), decreasing=T)) barplot(table(caixeta$local)) ## Cleveland dotplot par(mfrow=c(1,2)) plot(caixeta$h, ylab="Altura (m)", xlab="Observações") dotchart(caixeta$h, ylab="Observações", xlab="Altura (m)") par(mfrow=c(1,1)) ## Cleveland dotplot por grupo caixeta$localindex <- caixeta$local #objeto com indices para cores caixeta$localindex[caixeta$localindex=="chauas"] <- 1 caixeta$localindex[caixeta$localindex=="jureia"] <- 2 caixeta$localindex[caixeta$localindex=="retiro"] <- 3 caixeta$localindex <- as.numeric(caixeta$localindex) dotchart(caixeta$h, groups=factor(caixeta$local), col=caixeta$localindex) ## Identificando o outlier plot(caixeta$h) identify(caixeta$h) caixeta[362,] caixeta.2 <- caixeta[-362,] dim(caixeta) dim(caixeta.2) ## Boxplot boxplot(caixeta$h ~ caixeta$local, varwidth=T) #numero diferentes de observacoes em cada nivel ### Graficos univariados basicos ### ## Numa tela só boxplot, histograma, densidade e stripchart par(mfrow=c(2,2)) boxplot(aves.c$urubu) hist(aves.c$urubu) plot(density(aves.c$urubu)) stripchart(aves.c$urubu, method="stack") par(mfrow=c(1,1)) #boxplot representa 4-quantis, como em summary ou default do quantile quantile(aves.c$urubu) ## Stripchart ou dotplot stripchart(aves.c$urubu, method="stack") #para empilhar valores iguais ## Cleveland dotplot par(mfrow=c(2,1)) dotchart(aves.c$urubu, main="Urubu") boxplot(aves.c$urubu, horizontal=T) par(mfrow=c(1,1)) ## Histograma com os valores (funcao rug) hist(aves.c$urubu) rug(jitter(aves.c$urubu)) #representa obs hist(aves.c$urubu) # Sem o jitter rug(aves.c$urubu) hist(aves.c$urubu, breaks=20) ## Histograma com área = 1 e density kernel hist(aves.c$urubu, prob=T) lines(density(aves.c$urubu),col="blue", lwd=2) lines(density(aves.c$urubu, bw=5), col="red" ) # Diferentes tamanhos de banda lines(density(aves.c$urubu, bw=1.5), col="green" ) ## Adicionando uma curva da normal ## Primeiro: demonstração da funcao curve par(mfrow=c(1,2)) curve(expr=dnorm(x),from=-4, to=4, col="blue") # normal padronizada (mu=0, sigma=1) curve(expr=dnorm(x,mean=1),add=T, col="red") # aumentando mu curve(expr=dnorm(x),from=-4, to=4, col="blue") # normal padronizada de novo curve(expr=dnorm(x,sd=1.5),add=T, col="red") # aumentando sigma par(mfrow=c(1,1)) ## Sobrepondo a curva normal ao histograma hist(aves.c$urubu, prob=T) curve(expr = dnorm(x,mean=mean(aves.c$urubu),sd=sd(aves.c$urubu)), add=T, col="red") ## Linha de densidade no mesmo histograma lines(density(aves.c$urubu),col="blue") ## Sobrepondo a curva normal ao density plot(density(aves.c$urubu),col="blue", ylim=c(0,0.08)) curve(expr = dnorm(x,mean=mean(aves.c$urubu),sd=sd(aves.c$urubu)), add=T, col="red") ## Painel de histogramas require(lattice) histogram(~ caixeta$h | caixeta$local) densityplot(~ caixeta$h | caixeta$local) ## Grafico quantil-quantil ## Sorteio de 100 valores de uma normal com media = 30 e desvio-padrao = 3 ## Valores arredondados para 2 casa, e ordenados x <- sort(round(rnorm(100,30,3),2)) ## Inspecionando os 6 primeiros e ultimos valores head(x) tail(x) ## Media e dp mean(x) sd(x) ## Com a funcao qqnorm qqnorm(x) qqline(x, col="blue", lwd=2) ## Exemplo dados nao-normais qqnorm(aves.c$urubu) qqline(aves.c$urubu, col="red", lwd=2) ## Outro exemplo #ATENCAO: seu R precisa estar na versao mais recente para carregar a funcao qqPlot do pacote car require(car) qqPlot(caixeta$h) qqPlot(caixeta$h, distribution="pois", lambda=mean(caixeta$h)) qqPlot(caixeta$h, distribution="unif") qqPlot(caixeta$h, id.n=3) ############################################# ### Duas variaveis: tabulacao e agregacao ### ############################################# ## Tabelas ## ## Numero de fustes de cada especie por local table(caixeta$especie,caixeta$local) ## xtabs ## require(MASS) Titanic.df <- data.frame(Titanic) head(Titanic.df) ##Quanto sobreviventes por sexo? table(Titanic.df$Sex,Titanic.df$Survived) ## ops, não isso que queremos ## Precisamos da funcao xtabs xtabs(Freq~Sex+Survived, data=Titanic.df) prop.table(xtabs(Freq~Sex+Survived, data=Titanic.df), margin=1) xtabs(Freq~Class+Survived, data=Titanic.df) prop.table(xtabs(Freq~Class+Survived, data=Titanic.df), margin=1) ## tapply ## tapply(aves.c$carcara, aves.c$fisionomia, sum) tapply(aves.c$urubu, aves.c$fisionomia, sum) tapply(aves.c$urubu, aves.c$fisionomia, mean) ## aggregate para o dataframe caixeta ## ## Inspecionando o arquivo mais uma vez: head(caixeta, n=11) ## cada linha = um fuste (e nao uma arvore!) ## Criar um dataframe em que cada linha é uma arvore, com sua ## area basal total (soma da area de todo os fustes) ## calculando a area basal de cada fuste e incluindo no data frame caixeta caixeta$ab <- caixeta$cap^2/4*pi ## e agora criamos a planilha, com aggregate, ## somando as areas basais dos fustes caixeta.2 <- aggregate(caixeta$ab, by=list(local=caixeta$local, parcela=caixeta$parcela, arvore=caixeta$arvore, especie=caixeta$especie), FUN=sum) ## A nova variavel sai com nome "x" dim(caixeta.2) head(caixeta.2) ## mude o nome para ab ab <- names(caixeta.2)[5] str(caixeta.2) ## O mesmo para criar um data frame com a altura media por arvore temp <- aggregate(caixeta$h, by=list(local=caixeta$local, parcela=caixeta$parcela, arvore=caixeta$arvore, especie=caixeta$especie), FUN=mean) ## Inspecionando o objeto resultante head(temp) ## Antes de juntar a coluna de altura media, vamos nos certificar ## que ha correpondencia entre as linhas dos dois data frames sum(caixeta.2[,1:4]!=temp[,1:4]) # isto mostra que ha, entende pq? ## Adicione a coluna de altura media ao objeto caixeta.2 caixeta.2$h <- temp$x head(caixeta.2) ################################ ### Duas variaveis: graficos ### ################################ ## Lendo e conferindo a tabela esaligna <- read.table("esaligna.csv", header=T, sep=",", dec=".", as.is=T) str(esaligna) ## Grafico de dispersao plot(esaligna$ht~esaligna$dap) # Gráfico de dispersão com linha de suavização scatter.smooth(esaligna$ht~esaligna$dap, span=1/2) ## Scatterplot com boxplot do pacote car scatterplot(seriema~urubu, data= aves.c) ## Uso do argumento subset par(mfrow=c(1,3)) plot(seriema~urubu, data=aves.c, subset=fisionomia=="Ce", main="Ce") plot(seriema~urubu, data=aves.c, subset=fisionomia=="CC", main="CC") plot(seriema~urubu, data=aves.c, subset=fisionomia!="CL", main="CL") par(mfrow=c(1,1)) ## Exemplos da notacao de formula estatistica boxplot(urubu~fisionomia, data=aves.c) ############################ ### Quarteto de Anscombe ### ############################ data(anscombe)#carrega para a area de trabalho ls() #agora o objeto está no workspace names(anscombe) head(anscombe) apply(anscombe[1:4], MARGIN=2, FUN=mean) apply(anscombe[5:8], 2, mean) apply(anscombe[1:4], 2, var) apply(anscombe[5:8], 2, var) with(anscombe,cor(x1,y1)) with(anscombe,cor(x2,y2)) with(anscombe,cor(x3,y3)) with(anscombe,cor(x4,y4)) ## Graficos par(mfrow=c(2,2)) # 4 graficos em uma janela plot(y1~x1, data=anscombe) plot(y2~x2, data=anscombe) plot(y3~x3, data=anscombe) plot(y4~x4, data=anscombe) par(mfrow=c(1,1)) ## Linhas de regressao linear ## A fórmula da regressao é y= 3 + 0.5x para todos par(mfrow=c(2,2)) # 4 graficos em uma janela plot(y1~x1, data=anscombe, xlim=range(anscombe[,1:4]),ylim=range(anscombe[,5:8])) abline(lm(y1~x1, data=anscombe)) plot(y2~x2, data=anscombe, xlim=range(anscombe[,1:4]),ylim=range(anscombe[,5:8])) abline(lm(y2~x2, data=anscombe)) plot(y3~x3, data=anscombe, xlim=range(anscombe[,1:4]),ylim=range(anscombe[,5:8])) abline(lm(y3~x3, data=anscombe)) plot(y4~x4, data=anscombe, xlim=range(anscombe[,1:4]),ylim=range(anscombe[,5:8])) abline(lm(y4~x4, data=anscombe)) par(mfrow=c(1,1)) ##################### ### Multivariadas ### ##################### # Graficos condicionados do pacote car scatterplot(seriema~urubu|fisionomia, data= aves.c, lwd=2) ## Gráfico condicionado com o lattice ## library(lattice) xyplot(seriema~urubu|fisionomia, data= aves.c) coplot(seriema~urubu|fisionomia, data= aves.c) # Compare com a função padrão ## Outro tipo de gráfico condicionado xyplot(caixeta$h ~ caixeta$cap | caixeta$local) coplot(caixeta$h ~ caixeta$cap | caixeta$local) ## Tabelas para combinacoes de mais de duas variaveis ## with(caixeta,table(especie,parcela,local))# note uso da funcao with xtabs(Freq~Class+Survived+Sex, data=Titanic.df) prop.table(xtabs(Freq~Class+Survived+Sex, data=Titanic.df), margin=c(1,3)) ## Matrizes de correlacao ## cor(esaligna[,4:7]) cor(esaligna[,4:7], method="spearman") # coeficiente de Spearman ## Matriz de scatterplots ## pairs(esaligna[,4:7]) ## Matriz de similaridade ## aves.cf <- aggregate( aves.c[,2:4], by=list(fisio=aves.c$fisionomia), sum ) rownames(aves.cf) <- aves.cf$fisio aves.cf aves.cf.e <- dist(aves.cf[,2:4]) aves.cf.e plot(hclust(aves.cf.e), ) ## Exemplo de MDS com vegan ## library(vegan) plot(metaMDS(aves.cf[,2:4]),type="t") teste <- metaMDS(aves.cf[,2:4]) teste$species teste$points