<WRAP tabs>
  * [[02_tutoriais:tutorial4:start|Tutorial]]
  * [[01_curso_atual:exercicios4| Exercícios]]
  * [[03_apostila:05-exploratoria| Apostila]] 
</WRAP>
====== 4. Tutoriais de Análise Exploratória de Dados ======

Antes de iniciar uma análise propriamente dita, precisamos conhecer os dados e avaliar problemas com relação às etapas anteriores da pesquisa (coleta, processamento, estruturação e digitação dos dados brutos). As variáveis coletadas são aleatórias, ou seja, há diferentes fontes de variabilidade associadas a elas. Precisamos avaliar essa variabilidade e suas relações com as outras variáveis para poder tomar decisões embasadas sobre as análises que pretendemos realizar.

Dentre os principais **objetivos** de uma Análise Exploratória de Dados (AED) podemos listar os seguintes:

  * Detectar erros nos dados;
  * Detectar valores extremos (//outliers//);
  * Compreender a estrutura dos dados coletados;
  * Indicar a melhor distribuição para modelar a variabilidade nos dados;
  * Avaliar preliminarmente se os dados apresentam dependência/autocorrelação (espacial ou temporal);  
  * Avaliar se variáveis preditoras apresentam colinearidade;

Não é objetivo deste tutorial passar por todos esses tópicos e sim, apresentar algumas técnicas da linguagem R para auxiliá-lo a iniciar as análises exploratória dos dados. Essas técnicas passam por duas instrumentações básicas: análise numérica com estatísticas descritivas e análises gráficas, explorando a variabilidade das variáveis e sua relação com outras variáveis.
 

===== Conferindo Data Frames =====

Após a leitura dos dados no R é uma boa prática verificá-lo cuidadosamente e investigar as variáveis em busca de erros e incoerências. Nesse tutorial iremos explorar algumas das ferramentas básicas para explorar dados de variáveis.
<WRAP center round box 80%>

[[http://phdcomics.com/comics/archive.php?comicid=1612|{{ :02_tutoriais:tutorial4:phd01.gif?700 |}}]]

</WRAP>


Nesta sessão apresentamos duas situações bem comuns: variáveis de fatores com códigos errados e valores ''NA'' incorretos.

Vamos usar uma planilha com dados fictícios de um estudo de [[:dados::dados-aves|aves no cerrado]]. Baixe o arquivo do link no seu diretório de trabalho.

Vamos ler o arquivo de texto contendo os dados com o ''read.table''. Note que o padrão desse arquivo ''csv'' é o padrão da exportação do excel configurado em português.

<code rsplus>
aves <- read.table("aves_cerrado.csv", row.names = 1, header = TRUE, sep = ";", dec = ",", as.is = TRUE)
</code>

Logo após a leitura, garanta que os dados foram lidos corretamente, verificando a estrutura e outros atributos:

Verificação inicial do // data frame//
<code rsplus>
str(aves)
head(aves)
tail(aves)

</code>

Com isso já temos um diagnóstico de um possível problema: os ''NA''s são de fato observações faltantes?
Verificando as informações coletadas em campo descobrimos que esses campos são observações em que não houve avistamento da ave. Portanto, o valor correto é zero ((o problema inverso é comum e mais difícil de diagnosticar: colocar zero onde a observação não estava disponível!)).


==== Manipulando NA ====

Vamos indexar para verificar onde estão esses ''NA''s, usando a lógica básica de extração de registros em uma variável, no caso para ''urubu'': 

<code rsplus>
aves$urubu
aves$urubu == NA
</code>

O teste lógico para ''NA''s não funciona pois ele é uma palavra reservada e sua operação retorna ''NA''. Para o teste lógico específico de ''NA'' usamos a função ''is.na'', da seguinte forma:


<code rsplus>
is.na(aves$urubu)
which(is.na(aves$urubu))
</code>

O ''is.na'' retorna um vetor lógico que podemos usar para indexar o data frame ou apenas uma das suas variáveis:

<code rsplus>
aves[is.na(aves$urubu), ]
aves$fisionomia[is.na(aves$urubu)]
</code>

É possível também operar vários vetores lógicos para retornar os ''NA''s em mais de uma variável:

<code rsplus>
is.na(aves$urubu) | is.na(aves$carcara)
</code>

Acima o teste lógico com ''|'' entre vetores lógicos, opera nas posições utilizando a <wrap em>regra de equivalência</wrap>, retornando ''TRUE'' se pelo menos um dos vetores tiver, na mesma posição, o valor ''TRUE''. O teste retornará ''FALSE'' apenas se os vetores tiverem, na mesma posição, o valor ''FALSE''. Ao final, o vetor resultante terá ''TRUE'' quando ''urubu'' for ''NA'' <wrap em>ou</wrap> ''carcara'' for ''NA'' na posição equivalente. 

<code rsplus>
aves[is.na(aves$urubu) | is.na(aves$carcara) | is.na(aves$seriema), ]
</code>

Se indexarmos o objeto ''aves'' com a operação lógica ''is.na'' e ''|'' para os três vetores, ele retorna as linhas onde o teste é ''TRUE'', ou seja, onde existe algum ''NA'' entre as variáveis ''urubu'', ''caracara'' ou ''seriema''.
Vamos guardar esses registros em um objeto, atribuindo o resultado do código acima a um objeto:

<code rsplus>
avesNAs <- aves[is.na(aves$urubu) | is.na(aves$carcara) | is.na(aves$seriema), ]
</code>

Agora vamos sobrescrever os ''NA''s com o valor ''0'' nas posições associadas aos vetores lógicos produzidos pelo ''is.na''. Para garantir que está sobrescrevendo na posição onde está o erro, faça o processo em duas etapas, primeiro verifique o valor naquela posição antes de alterá-lo.

<code rsplus>
aves$urubu[is.na(aves$urubu)]
aves$urubu[is.na(aves$urubu)] <- 0
</code>

Fizemos para a variável ''urubu'' e podemos repetir o código para as outras variáveis. Uma maneira mais eficiente é usar o ''is.na'' em todo o objeto de uma única vez:


<code rsplus>
is.na(aves)
aves[is.na(aves)]
aves[is.na(aves)] <- 0
</code>

Para garantir que fizemos a substituição corretamente, vamos verificar, comparando as linhas que agora têm ''0''((Lembrando que antes não havia nenhum ''0'' nos dados, ou seja, as linhas que agora têm ''0'' são as linhas onde, anteriormente, sobrescrevemos algum elemento por ''0'')) com o objeto onde guardamos os registros contendo ''NA''s:

<code rsplus>
aves[aves$urubu == 0 | aves$carcara == 0 | aves$seriema == 0, ]
avesNAs
</code> 

==== Conferindo Fatores ====

As variáveis categóricas ou fatores podem ser conferidas com as funções ''table'' ou ''unique'' que retornam os níveis únicos, no primeiro caso fazendo a contagem e no segundo apenas os nomes únicos. Vamos aplicá-las para a variável ''fisionomia'':  


<code rsplus>
unique(aves$fisionomia)
table(aves$fisionomia)
</code>

O nível ''ce'' parece ser um erro de digitação já que todas as fisionomias deveriam ter 20 pontos amostrados. Para corrigir usamos a indexação usual de teste lógico:

<code rsplus>
aves$fisionomia == "ce"
aves[aves$fisionomia == "ce", ]
aves$fisionomia[aves$fisionomia == "ce"]
aves$fisionomia[aves$fisionomia == "ce"] <- "Ce"
table(aves$fisionomia)
</code>


Vamos agora converter esse vetor de caracteres para fator, ordenando as fisionomias da mais aberta para a mais fechada:

<code rsplus>
aves$fisionomia <- factor(aves$fisionomia, levels = c("CL", "CC", "Ce"))
</code>

Sempre é bom verificar novamente para ver se a estruturação e os valores do objeto estão corretos: 

<code rsplus>
str(aves)
</code>
Agora parece que tudo está correto!
===== Estatística Descritiva =====

Vamos usar o mesmo arquivo da sessão anterior para explorar as estatísticas descritivas básicas, começando pela média e pela mediana. Nós já usamos o ''apply'' para aplicar uma função a alguma dimensão de um objeto:


<code rsplus>
apply(X = aves[ ,2:4], MARGIN = 2, FUN = mean)
apply(aves[ ,c("urubu", "carcara", "seriema")], 2, median)
</code>

O que significaria a média ser muito diferente da mediana? A mediana é pouco sensível aos valores extremos e mais sensível à assimetria da distribuição do que a média. 

Podemos também fazer uma média truncada, usando o argumento ''trim'' da função ''mean''. Note que usamos uma outra indexação, com o mesmo resultado das anteriores, e que podemos incluir no ''apply'' os argumentos da função que iremos aplicar: 

<code rsplus>
apply(aves[ , -1], 2, mean, trim = 0.1)
</code>

O ''trim'' retira do cálculo da média os valores extremos com o corte definido pelo fator estipulado a partir das observações extremas. Nesse caso, retiramos 10% dos maiores valores e 10% dos menores. Como a média é muito sensível a valores extremos, se houver algum valor muito grande ou pequeno em relação ao resto, a média truncada seria bem diferente da média com todos os dados.
 
Os quantis também são uma forma de verificar se a distribuição dos valores é simétrica. O padrão da função ''quantile'' é retornar os quartis, que é a divisão dos dados em seus quartos depois de ordená-los: mínimo, 1/4, 1/2, 3/4 e máximo, sendo que a mediana é o segundo quartil onde os dados são divididos pela metade (menores e maiores). Quando usamos a função ''summary'' em um vetor de dados numéricos, esses valores também são apresentados:
<code rsplus>
quantile(aves$urubu) 
summary(aves$urubu)
</code>

O argumento ''probs'' da função ''quantile'' permite retornar o valor de qualquer quantil, vamos separar os dados em 10 quantis, com 10% dos dados cada um:

<code rsplus>
quantile(aves$urubu, probs = seq(from = 0, to = 1, by = 0.1))
</code>

A função ''summary'' ajuda a resumir essas estatísticas básicas para todas as variáveis numéricas ao mesmo tempo:
 
<code rsplus>
summary(aves[ ,-1])
</code>

===== Gráficos Univariados ======

Os gráficos são ferramentas importantes para avaliarmos as variáveis dos nossos dados. No tópico de gráficos vamos estudar mais a fundo as funções associadas a elaboração de gráficos no R, por enquanto vamos apenas aplicar algumas funções básicas sem nos preocuparmos muito com o acabamento dos gráficos. Normalmente a análise exploratória de dados é introspectiva, buscando entender as variáveis e não apresentá-las para alguma audiência.
Vejamos alguns gráficos básicos de diagnóstico de uma variável numérica, usando como exemplo o avistamento de urubus no Cerrado:

/*

<code rsplus>
caixeta <- read.table("http://ecor.ib.usp.br/lib/exe/fetch.php?media=dados:caixeta.csv", header = TRUE, sep = ",")

</code>
*/
==== Diagrama de Caixas ====

O diagrama de caixas (boxplot) é a representação gráfica dos quartis que a função ''quantile'' retorna.

<WRAP center round box 80%>
{{  :02_tutoriais:tutorial4:boxplot1-700x354.png?500  |  }}

</WRAP>

<code rsplus>
boxplot(aves$urubu)
</code>

==== Histograma e Gráfico de Barras ====

Os histogramas são usados para representar a frequência absoluta (contagem) ou densidade de variáveis contínuas e os gráficos de barras para variáveis categóricas ou discretas. Abaixo construímos as duas formas, mas antes repartimos nossa janela gráfica em quatro porções utilizando a função ''par'' que modifica os parâmetros gráficos e mudamos o parâmetro ''mfrow''.

<code rsplus>
par(mfrow = c(2, 2))
barplot(table(aves$urubu))
hist(aves$urubu)

</code>


O formato dos gráficos é muito diferente. Isso se deve ao fato do histograma agrupar as observações em intervalos, geralmente regulares. O número de intervalos vai influenciar muito o formato do gráfico, o R usa o algoritmo padrão da função ''nclass''. Vejam o que acontece se utilizarmos o argumento ''breaks'' do ''hist'': 

<code rsplus>
hist(aves$urub, breaks = 25)
</code>


Vamos usar a função ''stripchart'' para fazer mais um gráfico dessa variável e retornar o parâmetro ''mfrow'' para o seu padrão:


<code rsplus>
stripchart(aves$urubu, method="stack")
par(mfrow = c(1, 1))
</code>


No código abaixo incluímos os principais gráficos diagnósticos em uma única janela:


<code rsplus>
par(mfrow = c(2, 2))
boxplot(aves$urubu)
hist(aves$urubu)
barplot(table(aves$urubu))
stripchart(aves$urubu, method = "stack")
par(mfrow = c(1, 1))
</code>


==== Comparações com distribuições teóricas ====

Vamos agora retornar aos dados reais de [[dados:dados-caixeta|caixetais]] para explorar a variável numérica contínua ''cap'' das árvores:

<code rsplus>
caixeta <- read.table(file = "caixeta.csv", sep = ",", header = TRUE) #certifique-se que o arquivo esta em seu diretório de trabalho
summary(caixeta$cap)
par(mfrow = c(1, 2))
hist(caixeta$cap)
boxplot(caixeta$cap)
par(mfrow = c(1, 1))
</code> 

Essa distribuição de valores não parece nada com uma distribuição normal. Parece muito assimétrica com os dados concentrados nos valores menores. Além disso, tem uma longa cauda para os valores maiores. Várias características desses dados levam a uma distribuição como essa, as principais são:
  * o ''cap'' mínimo de inclusão foi 20 mm, ou seja não há como ter valores menores. Mesmo que pudesse, os valores estariam restritos a valores positivos e seria truncado no zero;
  * a estrutura de tamanho de populações estáveis de plantas tem a tendência a ter uma distribuição como essa de ''J'' invertido, devido às variações nas taxas vitais (crescimento, sobrevivência e reprodução) ao longo da sua ontogênese.

Uma forma de avaliar o acoplamento de um variável a uma distribuição normal é através dos gráficos ''qqnorm'' e da relação ''qqline'' que compara a distribuição dos valores dos quantis (pontos) da variável com os quantis de uma distribuição normal teórica (linha):

<code rsplus>
qqnorm(caixeta$cap)
qqline(caixeta$cap)
</code>  

Uma outra forma de fazer essa comparação é com o histograma de densidade, com o argumento ''freq = FALSE'' e utilizar a função ''curve'' que produz a curva da distribuição teórica:
 
<code rsplus>
hist(caixeta$cap, freq = FALSE)
curve(dnorm(x, mean = mean(caixeta$cap), sd = sd(caixeta$cap)), add = TRUE)
</code>  

Para resolver o problema de distribuição de valores assimétricos e positivos é muito comum utilizar a transformação ''log'', vamos explorar essa transformação:

<code rsplus>
par(mfrow = c(1, 2))
hist(caixeta$cap, freq = FALSE)
curve(dnorm(x, mean = mean(caixeta$cap), sd = sd(caixeta$cap)), add = TRUE)
hist(log(caixeta$cap), freq = FALSE)
curve(dnorm(x, mean = mean(log(caixeta$cap)), sd = sd(log(caixeta$cap))), add = TRUE)
par(mfrow = c(1, 1))
</code> 

Parece que a transformação logarítmica ajudou no acoplamento da variável à distribuição normal. 

<WRAP center round box 90%>

[[https://mathwithbaddrawings.com/2018/08/15/the-bubble-under-the-mathematical-rug/|{{  :02_tutoriais:tutorial4:normal.png?600  |}}]]

</WRAP>
 
/*
===== Exploração de uma Variável Categórica =====
Vamos usar um conjunto de dados de um inventário de árvores, que você pode baixar [[:dados:dados-caixeta| aqui]]. Leia com atenção a descrição deste conjunto de dados.

Vamos explorar a variável categórica nome da espécie, com a função ''table'':
 
<code rsplus>
caixeta <- read.csv("caixeta.csv", as.is = TRUE)
names(caixeta)
table(caixeta$especie)
</code>
Qual a diferença no resultado se usamos o comando abaixo?
<code rsplus>
sort(table(caixeta$especie), decreasing = TRUE)
</code>

O resultado da função ''table'' pode ser fornecido à função de gráficos de barras:
<code rsplus>
barplot(sort(table(caixeta$especie), decreasing = TRUE))
barplot(table(caixeta$local))
</code>


===== Variações do Histograma =====
Voltando ao objeto criado no tutorial [[:02_tutoriais:tutorial4:start#conferindo_data_frames|Conferindo Data Frames]], vamos fazer algumas variações de histogramas do número de avistamentos de urubus:
Você pode acrescentar marcas (traços) indicando a posição de cada observação no eixo x.
<code rsplus>
## Histograma com os valores (função rug)
hist(aves.c$urubu)
rug(jitter(aves.c$urubu))
</code>
O que acontece se você omite a função ''jitter'' neste caso? Por que?

Agora vamos fazer um histograma re-escalonado de modo que as áreas das barras somem um. Com isto, podemos sobrepor ao histograma um ajuste não paramétrico de densidade probabilística, que também mantém área um:

<code rsplus>
hist(aves.c$urubu, prob = TRUE)
lines( density(aves.c$urubu), col = "blue" )
</code>

Também sobre este histograma podemos sobrepor a curva normal. Para os parâmetros da normal, usamos a média e o desvio-padrão da amostra.
<code rsplus>
hist(aves.c$urubu, prob = TRUE)
curve(expr = dnorm(x, mean = mean(aves.c$urubu), sd = sd(aves.c$urubu)), add = TRUE, col ="red")
</code>

Por fim, vamos sobrepor a curva empírica de densidade probabilística com a curva normal:
<code rsplus>
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")
</code>

O que estes gráficos revelam sobre a distribuição do número de avistamentos de urubus neste estudo fictício?
*/
===== Análise exploratória Bivariada =====
 
A relação entre duas ou mais variáveis categóricas pode ser explorada com tabelas cruzadas, por exemplo:

<code rsplus>
table(caixeta$especie,caixeta$local)
</code>


Quando temos uma variável categórica (fator) e uma numérica, as funções ''aggregate'' e ''tapply'' são muito úteis. 

A função ''aggregate'' é o equivalente das tabelas dinâmicas das planilhas eletrônicas: ela aplica uma função a partes dos dados separadamente. Por exemplo, para obter um //data frame// com a altura média dos fustes de cada espécie de árvore por local você executa o comando:

<code rsplus>

caixetaH <- aggregate(caixeta$h, by = list(local = caixeta$local, especie = caixeta$especie), FUN=mean)

</code>

A função ''aggregate'' recebe uma ou mais variáveis contínuas e aplica a ela(s) a função que é recebida no argumento ''FUN'' agrupadas pelas variáveis que foram colocadas no argumento ''by'' organizadas em uma lista. Veja a estrutura do objeto resultante:

<code rsplus>
str(caixetaH) 
</code>

A função ''tapply'' faz a mesma operação, mas organiza a objeto resultante em uma ''matrix'', veja como ficam organizadas as médias de altura nesse caso:

<code rsplus>
caixAlt <- tapply(caixeta$h, INDEX = list(caixeta$especie,especie=caixeta$local), FUN = mean) 
str(caixAlt)
class(caixAlt)
</code>

<code rsplus>
dim(caixetaH)
dim(caixAlt)
</code>

Os dados originais continha 1027 observações, agora o ''caixetaH'' tem 48, pois compilou os dados de altura em sua média para cada espécie em cada uma das localidades. O objeto ''caixAlt'' por sua vez tem 43 linhas e 3 colunas que contém as mesmas médias. Como é possível?  \\
No primeiro caso nós temos 3 variáveis: ''local'', ''especie'' e ''x'' (média de altura). No ''caixAlt'' temos as linhas representando as ''espécies'', as colunas representando as ''localidades'' e os valores no interior da ''matrix'' como as médias das alturas.
 

Consulte a ajuda da função ''aggregate'' e experimente outras combinações de fatores e funções, com este conjunto de dados.

/*
===== xtabs =====
Crie um objeto com este [[:dados::dados-titanic|arquivo]] e faça as seguintes tabulações:
<code rsplus>
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)
</code>

Por que usamos a função ''xtabs'' neste caso e não a função ''table''? P.ex:
<code rsplus>
table(Titanic.df$Sex, Titanic.df$Survived)
</code>
*/
===== Gráficos Bivariados =====

O gráfico de dispersão nos permite avaliar a relação entre duas variáveis numéricas:

<code rsplus>
plot(x = caixeta$cap, caixeta$h)
</code>


A função ''pairs'' auxilia na avaliação da relação entre múltiplos pares de variáveis, em um painel:

<code rsplus>
str(aves)
pairs(aves[ , c("urubu", "carcara", "seriema")])
</code>


A notação de fórmula estatística do R, usa o símbolo ''~'' para indicar a relação entre duas variáveis. 

<code rsplus>
fun(y ~ x, data = dados)
</code>

  * ''fun'' = função aplicada no objeto ''dados''
  *  ''y'' = variável resposta
  *  ''x'' = variável preditora
  * ''dados'' = um data frame contendo os vetores x e y

Esta notação foi criada para os modelos estatísticos, como a regressão linear, mas foi estendida para várias funções gráficas no R onde as relações entre as variáveis pode ser representada. Essa notação facilita muito a produção de gráficos exploratórios da relação entre variáveis. Os comandos abaixo utilizam esta notação:


<code rsplus>
plot(h ~ cap, data = caixeta)
plot(h ~ cap, data = caixeta, subset = local == "jureia")
</code>

A notação em fórmula nos gráficos tornam o código mais sintético e permitem algumas formulações interessantes, por exemplo entre variáveis contínuas e categóricas:

<code rsplus>
boxplot(cap ~ local, data = caixeta)
</code>


Esta fórmula pode ainda incluir um fator condicionante, que aplica a relação proposta dentro de cada nível dos condicionais:


<code rsplus>
fun(y ~ x | z, data = dados)
</code>
  * ''z'' = variável condicionantes

O pacote ''lattice'' implementa esta ideia para gráficos:

<code rsplus>
library(lattice)
xyplot(h ~ cap | local, data = caixeta)
</code>


===== O quarteto de Anscombe =====

O pacote ''datasets'' tem vários conjuntos de dados que se tornaram clássicos da estatística. Vamos analisar um dos mais usados para treino de análises exploratórias. Como o R já carrega o ''datasets'' por padrão, podemos utilizar os dados diretamente ((a função ''data'' carrega os dados do pacote no nosso ''workspace'', mas isto não é necessário já que ele está disponível na sessão a partir do ''environment'' ''datasets'')):

<code rsplus>
str(anscombe)

</code>

Este objeto é composto de 4 pares de variáveis, nomeadas ''x1'' a ''x4'' (variáveis independentes ou preditoras) e ''y1'' a ''y4'' (variáveis dependentes ou resposta):
<code rsplus>
names(anscombe)
</code>

Compare as médias de cada uma das variáveis. Para isto, use a função ''apply'' para aplicar a função ''mean'' a cada coluna((argumento ''MARGIN=2'')) deste data frame:
<code rsplus>
apply(anscombe[1:4], MARGIN = 2, FUN = mean)
apply(anscombe[5:8], 2, mean)
</code>

Faça o mesmo para obter as variâncias:

<code rsplus>
apply(anscombe[1:4], 2, var)
apply(anscombe[5:8], 2, var)
</code>

A pergunta principal para este conjunto de dados é se há relação entre cada variável ''x'' e ''y''. Isso pode ser avaliado com o [[http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient|coeficiente de correlação de Pearson]], que vai de zero (nenhuma correlação) a um (positivo ou negativo, correlação perfeita).

<code rsplus>
cor(anscombe$x1, anscombe$y1)
</code>

O código acima pode ser aplicado através da função ''with'' que simplifica a digitação das variáveis assim como o argumento ''data'' na função ''plot'':

<code rsplus>
with(anscombe, cor(x1, y1))
with(anscombe, cor(x2, y2))
with(anscombe, cor(x3, y3))
with(anscombe, cor(x4, y4))
</code>

Uma outra forma de avaliar a relação é através dos coeficientes da reta de um modelo linear entre as duas variáveis. Essa relação é obtida pela função ''coef'' aplicada ao um objeto de modelo linear ((isso será alvo das aulas de modelo linear)). Vamos avaliar esse coeficientes para os quatro conjuntos de variáveis:

<code rsplus>
(coef1 <- coef(lm(y1 ~ x1, data = anscombe)))
(coef2 <- coef(lm(y2 ~ x2, data = anscombe)))
(coef3 <- coef(lm(y3 ~ x3, data = anscombe)))
(coef4 <- coef(lm(y4 ~ x4, data = anscombe)))
</code> 

Perceba que todos os coeficientes estão muito próximos a ''3.00'' e ''0.50''.
Vamos fazer os gráficos de dispersão entre as variáveis ''x'' e ''y'' para os quatro casos e vamos incluir a relação linear para cada par com a função ''abline'' que grafa essa relação. Nossa única inclusão no ''plot'' foram os argumentos ''xlim'' e ''ylim'' para que todos os gráficos tenham as mesmas escalas nos eixos ((veja a documentação do ''plot'')):


/*
<code rsplus>
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))
</code>
*/

<code rsplus>
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(coef1)
plot(y2 ~ x2, data = anscombe,
     xlim = range(anscombe[,1:4]), ylim = range(anscombe[,5:8]))
abline(coef2)
plot(y3 ~ x3, data = anscombe,
     xlim = range(anscombe[,1:4]), ylim = range(anscombe[,5:8]))
abline(coef3)
plot(y4 ~ x4, data = anscombe,
     xlim = range(anscombe[,1:4]), ylim = range(anscombe[,5:8]))
abline(coef4)
par(mfrow=c(1,1))
</code>

Este conjunto de dados foi criado pelo estatístico Frank Anscombe para demonstrar a importância da análise visual de dados (veja mais detalhes [[http://en.wikipedia.org/wiki/Anscombe's_quartet|aqui]]). Em nossas análises exploratórias não gráficas, vimos que as médias, desvios, correlações e coeficientes de reta são praticamente idênticos, o que indicaria que os dados são similares. Note como são diferentes, apesar das estatísticas similares: o primeiro é uma relação que parece estar adequada às premissas dos modelos lineares, o segundo mostra uma clara relação não linear entre as variáveis, já o terceiro tem um dado influente que promove uma inclinação que não acompanha o conjunto dos dados e promove além da mudança da relação uma não homogeneidade da variância, por fim, o quarto mostra também um ponto muito influente e com alta alavancagem que define a relação entre as variáveis e, caso seja retirado da amostra, as variáveis ''y4'' em função de ''x4'' não apresentam nenhuma relação! 

===== Leia mais sobre análise exploratória de dados =====


Artigo com um protocolo de análise exploratória de dados:
  *[[http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/pdf]]

Capítulo do livro "Analysing Ecological Data" sobre análise exploratória de dados:
  *{{:02_tutoriais:tutorial4:zuur_cap_4_exploration_statistics.pdf|}}