Gráficos de barras (e um pouco de ANOVA) em R

(Este post pressupõe um conhecimento básico de R. Tenho um post introdutório, caso queira, aqui).

Às vezes a melhor forma de apresentar os nossos resultados é um gráfico de barras, por exemplo quando queremos comparar entre diferentes níveis de alguns tratamentos. Por exemplo, em um artigo [recém-aceito, uhul!], comparamos o tempo de sobrevivência de ninhos artificiais entre duas localidades (borda e interior) e duas alturas (alto e baixo). O gráfico resultante foi o seguinte:

fig1.png

Nele vemos simultaneamente a relação do tempo de sobrevivência com o local onde o ninho está e com a sua altura. Resumidamente, apenas os ninhos baixos localizados no interior têm uma taxa de sobrevivência maior. As barras são intervalos de confiança, calculados por bootstrap.

Mas, pra deixar esse post mais generalizável e porque agora a Elsevier detêm os direitos autorais sobre o artigo, melhor não arriscar, vamos imaginar que estamos trabalhando com… hummm… Comprimento de cauda de sereias de mangue, obviamente depois de passar pelos comitês de ética para trabalhos com seres humanos, animais não-humanos, e seres não-humanos fantásticos. Digamos que foram escolhidas cinco áreas em um manguezal, cada área com três tipos de vegetação ou estratos (vegetação recente, intermediária e antiga), e quinze sereias de cada estrato de cada área tiveram suas caudas medidas, por pesquisadoras responsáveis e evitando causar qualquer tipo de desconforto. Todas as sereias participaram do trabalho de forma voluntária.

Não tendo os dados reais deste estudo, vamos simular um resultado possível. Digamos que não há difereça entre os tipos de vegetação, mas há diferença entre as áreas amostradas.

Vamos primeiro criar as nossas variáveis explanatórias:

area = rep(c("Area1", "Area2", "Area3", "Area4", "Area5"), each=45)
estrato = rep(rep(c("Recente", "Intermediaria", "Antiga"), each=15), times=5)
dados = data.frame(area, estrato)

E agora vamos simular os valores. Vamos considerar um comprimento médio de 80 cm para a área 1, 90 cm para a área 2 etc, com desvio-padrão de 10 cm e distribuição normal.

Primeiro, vamos definir o valor médio correspondente a cada sereia, com base na área onde ela vive:

area.num = as.numeric(as.factor(area))
area.comp = 70 + area.num*10

A partir deste vetor de médias, geramos os de comprimento. Temos um total de 5 areas * 3 estratos * 15 sereias = 225 valores de comprimento de cauda.

set.seed(17)
cauda = rnorm(n=225, mean=area.comp, sd=10)

Usei o comando set.seed para que os valores sejam sempre os mesmos nas simulações.

Uma forma de visualizar estes valores individualmente é fazendo um simples gráfico de dispersão:

plot(cauda, xlab="Ordem das observações", ylab="Comprimento da cauda (cm)", col=area.num, pch=14+area.num)

fig2.png

Neste gráfico, o eixo X não representa nada, sendo apenas a ordem em que as observações estão na planilha; o eixo Y é o comprimento de cauda; e as diferentes cores e formas dos símbolos representam as cinco áreas.

Agora, juntamos estes valores ao nosso objeto de dados:

dados$cauda = cauda

Se quisermos olhar a estrutura deste objeto, usamos o comando str(dados), que fornece o seguinte resultado:

str(dados)
'data.frame': 225 obs. of 3 variables:
$ area : Factor w/ 5 levels "Area1","Area2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ estrato: Factor w/ 3 levels "Antiga","Intermediaria",..: 3 3 3 3 3 3 3 3 3 3 ...
$ cauda : num 69.8 79.2 77.7 71.8 87.7 ...

Ou seja, temos duas variáveis categóricas, ou fatores, com 5 e 3 níveis respectivamente, e uma variável numérica.

Fazendo o gráfico

O R faz automaticamente boxplots mostrando a relação do tamanho de cauda com a área e o estrato:

par(mfrow=c(1,2))
plot(cauda~area+estrato, data=dados, ylab="Comprimento da cauda (cm)")

fig3.png

O argumento ylab nos dá o nome do eixo Y.

Esse gráfico já nos mostra o que queremos saber. Eu faria outros gráficos, mas vamos direto pro gráfico de barras, afinal, este é o título do post!

(Provavelmente há formas mais fáceis de fazê-lo. Mas, sou preguiçoso e não fui atrás de descobrí-las, rs.)

(Há um texto talvez melhor que esse, em inglês, aqui, e provavelmente outros, também melhores que esse, em português. Já falei que sou preguiçoso? 🙂 )

Em R, para o gráfico de barras precisamos fornecer os que corresponderão às barras. Ou seja, não podemos entrar direto com os dados, precisamos entrar com os valores médios e, se quisermos, com os valores de erros também.

Para isso, podemos usar a função aggregate, que aplica alguma função de interesse para cada nível das nossas variáveis categóricas:

caudas.mean = aggregate(cauda~area+estrato, data=dados, FUN=mean)
caudas.sd = aggregate(cauda~area+estrato, data=dados, FUN=sd)

O objeto caudas.mean nos mostra o valor médio em cada combinação de área e estrato, e o caudas.sd o desvio padrão:

caudas.mean
area estrato cauda
1 Area1 Antiga 74.80873
2 Area2 Antiga 88.92437
3 Area3 Antiga 96.64252
4 Area4 Antiga 110.63772
5 Area5 Antiga 122.51882
6 Area1 Intermediaria 82.50872
7 Area2 Intermediaria 86.00810
8 Area3 Intermediaria 103.61631
9 Area4 Intermediaria 109.04462
10 Area5 Intermediaria 122.40882
11 Area1 Recente 84.44749
12 Area2 Recente 94.20509
13 Area3 Recente 98.33244
14 Area4 Recente 109.44055
15 Area5 Recente 115.26678

caudas.sd
area estrato cauda
1 Area1 Antiga 12.201999
2 Area2 Antiga 11.271872
3 Area3 Antiga 13.297709
4 Area4 Antiga 8.336021
5 Area5 Antiga 9.391670
6 Area1 Intermediaria 7.204655
7 Area2 Intermediaria 11.431133
8 Area3 Intermediaria 10.402698
9 Area4 Intermediaria 13.790228
10 Area5 Intermediaria 7.771806
11 Area1 Recente 8.263902
12 Area2 Recente 11.377898
13 Area3 Recente 8.867977
14 Area4 Recente 10.581011
15 Area5 Recente 12.195718

Agora começa a parte mais chatinha.

Queremos um gráfico com cinco conjuntos de três barras. Os conjuntos são as áreas, as barras são os estratos. Para isso, precisamos transformar estes nossos objetos em matrizes.

caudas.mean.bp = matrix(caudas.mean$cauda, nrow=3, ncol=5, byrow=T)

Esse comando coloca os valores de comprimento de cauda em ordem em uma matriz; o argumento byrow=T faz com que eles sejam inseridos por linhas, não por colunas:

caudas.mean.bp
[,1] [,2] [,3] [,4] [,5]
[1,] 74.80873 88.92437 96.64252 110.6377 122.5188
[2,] 82.50872 86.00810 103.61631 109.0446 122.4088
[3,] 84.44749 94.20509 98.33244 109.4406 115.2668

Vamos agora dar nomes às linhas e colunas:

colnames(caudas.mean.bp) = c("Area1", "Area2", "Area3", "Area4", "Area5")
rownames(caudas.mean.bp) = c("Antiga", "Intermediaria", "Recente")

O resultado:

caudas.mean.bp
Area1 Area2 Area3 Area4 Area5
Antiga 74.80873 88.92437 96.64252 110.6377 122.5188
Intermediaria 82.50872 86.00810 103.61631 109.0446 122.4088
Recente 84.44749 94.20509 98.33244 109.4406 115.2668

Façamos o mesmo para os desvios:

caudas.sd.bp = matrix(caudas.sd$cauda, nrow=3, ncol=5, byrow=T)
colnames(caudas.sd.bp) = c("Area1", "Area2", "Area3", "Area4", "Area5")
rownames(caudas.sd.bp) = c("Antiga", "Intermediaria", "Recente")

E agora pra fazer o gráfico:

barplot(caudas.mean.bp, beside=T, legend.text=T, args.legend=list(x=7), ylab="Comprimento de cauda (cm)")

Na função acima, o argumento beside=T diz que as colunas ficarão uma do lado da outra, e não empilhadas; legend.text=T diz que teremos uma legenda explicando as cores; args.legend=list(x=7)) especifica a localização da legenda no eixo X. Cheguei ao valor 7 por tentativa e erro.

E olhem só que coisa bonita!
fig4.png

Mas e os erros?

Inserir os erros é um pouco mais trabalhoso. Precisamos especificar exatamente de onde as linhas partem e onde elas terminam, em termos de coordenadas X e Y.

As barras de erro partem dos nossos valores médios, ou seja, os valores do objeto caudas.mean.bp. E terminam nestes valores ± os desvios-padrões do nosso objeto caudas.sd.bp

As coordenadas X são os centros das nossas barras, produzidos pelo barplot. E portanto precisamos extraí-las do próprio barplot. Para isso, criamos um objeto:

caudas.X = barplot(caudas.mean.bp, beside=T)

Podemos fechar ele se quisermos:

dev.off()

O nosso objeto caudas.X tem as coordenadas X, com a mesma estrutura que as nossas matrizes:

caudas.X
[,1] [,2] [,3] [,4] [,5]
[1,] 1.5 5.5 9.5 13.5 17.5
[2,] 2.5 6.5 10.5 14.5 18.5
[3,] 3.5 7.5 11.5 15.5 19.5

Então fazemos novamente o nosso barplot, adicionando segmentos, ou setas, para os erros:

barplot(caudas.mean.bp, beside=T, legend.text=T, args.legend=list(x=7), ylab="Comprimento de cauda (cm)")
arrows(x0=caudas.X, x1=caudas.X, y0=caudas.mean.bp-caudas.sd.bp, y1=caudas.mean.bp+caudas.sd.bp, code=3, angle=90, length=0.1)

fig5.png
O argumento code=3 diz pra colocar uma seta nos dois fins do segmento, e o angle=90 diz que é pros lados da seta serem perpendiculares à linha (como uma barra de erro deve ser). x0, x1, y0 e y1 dizem as coordenadas X e Y iniciais e finais de cada segmento. Cada segmento vai de (média – desvio) até (média + desvio), com as coordenadas X dadas pelo próprio barplot. O argumento length=0.1 é pras barras horizontais ficarem menores; também cheguei no valor de 0.1 por tentativa e erro.

Parece um tanto complicado, né? Também acho! Não conheço uma forma mais simples de fazer estes gráficos em R, ser ser usando algum pacote como ggplot ou lattice (os quais não sei usar). Mas, cá está. O esquema geral é esse.

E a significância?

Agora que fizemos o gráfico, podemos querer também ver se há diferenças significativas… Dá pra usar uma ANOVA bifatorial. Isso é simples:

caudas.aov = aov(cauda~area*estrato, data=dados)

Ou seja, estamos testando se há um efeito da área, do estrato, e da interação entre as duas variáveis sobre o comprimento de cauda. Para testar se as premissas da análise não foram violadas, fazemos gráficos de validação:

par(mfrow=c(2,2))
plot(caudas.aov)

fig6.png
Não há padrões residuais, os resíduos se ajustam a uma distribuição normal, e não há outliers ou pontos influentes muito chamativos. Então parece que esta tudo em ordem. Ainda bem, afinal os dados foram simulados a partir de uma distribuição normal!

O resultado, obtido com o comando summary(caudas.aov):

Df Sum Sq Mean Sq F value Pr(>F)
area 4 44085 11021 97.885 <2e-16 ***
estrato 2 171 86 0.761 0.4687
area:estrato 8 2062 258 2.289 0.0227 *
Residuals 210 23645 113
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Há um efeito significativo da área e, estranhamente, também uma interação entre área e estrato. Como simulei o mesmo efeito de área pra todos os estratos, essa interação significativa não é esperada. É possível que seja um erro do tipo I (erroneamente rejeitar uma hipótese nula verdadeira).

Mas beleza, sabemos que há um efeito da área. Mas quais áreas diferem entre si? Para isso, fazemos um teste de Tukey sobre o resultado da ANOVA:

TukeyHSD(caudas.aov)

Ele dá um resultado bem longo, mostrando 1) quais áreas diferem de quais áreas, 2) quais estratos diferem de quais estratos, 3) quais combinações de área e estrato diferem de quais outras combinações. Resumidamente, todas as áreas diferem de todas as outras áreas (p < 0.0007), nenhum estrato difere de nenhum outro estrato, e vou deixar a interpretação das interações como lição de casa porque estou com preguiça de intepretar porque acho mais didático. 🙂

Espero que tenham se divertido! Um dia escrevo mais sobre como interpretar essas tabelas de ANOVA e teste de Tukey, ou não.

Anúncios

4 pensamentos sobre “Gráficos de barras (e um pouco de ANOVA) em R

  1. Pingback: Ensinando R… – Mais Um Blog de Ecologia e Estatística

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s