## Scripts da aula de analise exploratoria de dados

################################################################################
## Primeiro exemplo: uso da funcão is.na e substituicão de valores
################################################################################
## 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)
## primeiro problema: os NAs sao de fato observacões faltantes?
## Verificando a planilha descobrimos que na verdade sao pontos em que nao houve avistamento.
## Portanto, o valor correto e zero
## Quais sao os registros com este problema?
## Vamos ver para os urubus
aves.c[aves.c$urubu==NA,] ## Ops, nao funciona
## O teste logico Para NA e feito pela funcao is.na:
is.na(aves.c)
is.na(aves.c$urubu)
## Retomando entao:
aves.c[is.na(aves.c$urubu)==T,]
## E para NA em uma das tres aves
aves.c[is.na(aves.c$urubu)==T|is.na(aves.c$carcara)==T|is.na(aves.c$seriema)==T,]
## Vamos guardar este pedaco da planilha num objeto temporario
temp1 <- aves.c[is.na(aves.c$urubu)==T|is.na(aves.c$carcara)==T|is.na(aves.c$seriema)==T,]
## E agora vamos corrigir estes valores
aves.c$urubu[is.na(aves.c$urubu)==T] <- 0
## Que e o mesmo que
aves.c[is.na(aves.c$urubu)==T,2] <- 0
##Ou ainda
aves.c[is.na(aves.c[,2])==T, 2] <- 0
## Continuando, para as outras aves
aves.c$carcara[is.na(aves.c$carcara)==T] <- 0
aves.c$seriema[is.na(aves.c$seriema)==T] <- 0
## Verificando
aves.c[is.na(aves.c$urubu)==T|is.na(aves.c$carcara)==T|is.na(aves.c$seriema)==T,]
## 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!
##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)
## 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)

################################################################################
## Resumo estatistico: medias, media truncada e mediana
################################################################################
summary(aves.c[,2:4])
sd(aves.c[,2:4])
apply(aves.c[,2:4],2,median)
mean(aves.c[,2:4])
mean(aves.c[,2:4], trim=0.1)
## Quantis
quantile(aves.c$urubu) ## O mesmo que o retornado pelo summary, veja:
summary(aves.c$urubu)
## Mas vc pode mudar os quantis
quantile(aves.c$urubu, probs= seq(from=0,to=1,by=0.1))

################################################################################
## Tabelas de variaveis categoricas
################################################################################
caixeta <- read.csv("caixeta.csv", as.is=T) ##arquivo caixeta.csv deve estar no diretorio de trabalho
names(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))

################################################################################
## Graficos univariados
################################################################################

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

## Histograma com diferentes larguras de barras
par(mfrow=c(2,2))
hist(aves.c$urubu, main="Default")
hist(aves.c$urubu, breaks=seq(0,max(aves.c$urubu),length=5), main="Cinco break-points")
hist(aves.c$urubu, breaks=seq(0,max(aves.c$urubu),length=8), main="Oito break-points" )
hist(aves.c$urubu, breaks=seq(0,max(aves.c$urubu),length=10), main="Dez break-points")
par(mfrow=c(1,1))

## Algumas variacoes sobre o histograma
## Histograma com os valores (funcao rug)
hist(aves.c$urubu)
rug(jitter(aves.c$urubu))
## Histograma com área = 1 e density kernel
hist(aves.c$urubu, prob=T)
lines(density(aves.c$urubu),col="blue")
## Adicionando uma curva da normal aos graficos
## Usamos a funcao curve, que traca a curva de uma f(x), num intervlo especificado de x
## Demonstracao do curve
## Para a funcao da normal padronizada (media zero e desvio-padrao =1) no intervalo -2.5 a 2.5: 
curve(expr=dnorm(x),from=-2.5, to=2.5)
## Agora adionamos no mesmo grafico a funcao de densidade probabilistica de t
curve(expr=dt(x, df=5), add=T, col="red", lty=2) ## note argumento add=T
## Usando curve para sobrepor 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")
## 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")

################################################################################
## Exemplo para o qqplot
################################################################################
##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 5 primeiros e ultimos valores
x[1:5]
x[95:100]
## Calculo do percentil de cada valor
px <- order(x)/100
## Media e dp
mean(x)
sd(x)
## Quantis correspondentes aos percentis, para uma normal com media e dp tomadas da
## media e dp da amostra
q.norm.x <- qnorm(px,mean=mean(x),sd=sd(x))
## Juntando valores, percentis e quantis previstos em um dataframe, para facilitar a visualizacao
qq.plot.x <- data.frame(x=x, percentil=px, q.norm=q.norm.x)
qq.plot.x[1:5,]
qq.plot.x[95:100,]

## Grafico
plot(x~q.norm, data=qq.plot.x, xlab="Quantis Esperados", ylab="Valores Observados")
abline(0,1, col="red") # relacao esperada, caso os dados venhamd e uma populacao normal

## A funcao qqnorm ja faz isto de uma vez para voce:
qqnorm(x)
qqline(x)
qqnorm(aves.c$urubu)
qqline(aves.c$urubu)


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

##SEGUNDA PARTE: RELACOES ENTRE VARIAVEIS

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

## Tabulacao cruzada de dois ou mais fatores
## Numero de fustes de cada especie por local
table(caixeta$especie,caixeta$local)

##xtabs: tabulacao de dados de frequencia
## Dataframe dos sobreviventes dos sobreviventes e mortos do Titanic
head(Titanic.df)
##Quanto sobreviventes por sexo?
table(Titanic.df$Sex,Titanic.df$Survived) ## ops, nao e 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)
## E para combinacoes de mais de duas variaveis
xtabs(Freq~Class+Survived+Sex, data=Titanic.df)
prop.table(xtabs(Freq~Class+Survived+Sex, data=Titanic.df), margin=c(1,3))
      

## tapply: resumo de uma variavel numerica, separada por niveis de um ou mais fatores
tapply(aves.c$carcara,aves.c$fisionomia, sum)
tapply(aves.c$urubu,aves.c$fisionomia, sum)
tapply(aves.c$urubu,aves.c$fisionomia, mean)


## "Tabelas dinamicas": funcao aggregate
## aggregate para o dataframe caixeta
## Criar data.frame com altura media dos fustes por especie e por local
##para mostrar o aruivo mais uma vez
names(caixeta)
caixeta.alt <- aggregate(caixeta$h, by=list(local=caixeta$local,especie=caixeta$especie),
                         FUN=max)
head(caixeta.alt)
## Criar um dataframe em que cada linha e uma arvore, com sua area basal total
##(soma da area de todo os fustes)
## calculando a area basal de cada fuste
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)

## Exemplos de Graficos bivariados

## boxplot
boxplot(urubu~fisionomia, data=aves.c)

##espalhagrama
esalig <- read.csv("esaligna.csv")
plot(ht~dap, data=esalig)
plot(ht~dap, data=esalig, subset=classe=="a")
scatter.smooth(esalig$ht~esalig$dap)
scatter.smooth(esalig$ht~esalig$dap, span=1/2)

##Graficos com parte dos dados: argumento subset
## Um grafico para cada fisionomia
plot(seriema~urubu, data=aves.c, subset=fisionomia=="Ce")
plot(seriema~urubu, data=aves.c, subset=fisionomia=="CC")
plot(seriema~urubu, data=aves.c, subset=fisionomia!="CL")


##pairs
## Matriz de espalhagramas das medidas das arvores no dataframe esalig
pairs(esalig[,4:8])
## Que e a expressao grafica da matriz de correlacoes:
round(cor(esalig[,4:8]), 2)


## Graficos condicionados com o pacotec lattice
library(lattice) # carregue o pacote (basta fazer isto uma vez na secao)
xyplot(folha~tronco|classe, data=esalig)
## um data.frame com as duas especies mais abundantes
caixeta.abund <- caixeta[caixeta$especie=="Tabebuia cassinoides"|caixeta$especie=="Myrcia sulfiflora",]
bwplot(h~local|especie, data=caixeta.abund)
xyplot(h~cap|especie+local, data=caixeta.abund)
xyplot(seriema~urubu|fisionomia, data= aves.c)