Curso de R no IPÊ
Script do dia 3.1: ANOVA
Marcos V. C. Vital e participantes da turma de 2014
Para mais material, acesse http://marcosvital.wordpress.com/
Este documento é o terceiro de uma série que nasceu durante as atividades do “Curso de introdução ao uso do software R e suas aplicações estatísticas nas ciências biológicas”, realizado no IPÊ em dezembro de 2014.
Neste terceiro guia nós vamos: realizar uma ANOVA, brincar com mais algumas funcionalidades gráficas do R e, de quebra, dar uma olhadinha em como se faz um teste de Kruskal Wallis.
Como sempre, todo o documento foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais, acesse:
Antes de começar, carregue os dados necessários:
setwd("C:/R/Curso IPÊ")
dados<-read.table("cultivo.txt", h=T)
attach(dados)
summary(dados)## UA pragas cultivo
## Min. : 1.00 Min. :12.00 agroflorestal:14
## 1st Qu.:15.25 1st Qu.:38.00 ausente :15
## Median :29.50 Median :46.00 mata :15
## Mean :29.50 Mean :48.21 refugio :14
## 3rd Qu.:43.75 3rd Qu.:58.00
## Max. :58.00 Max. :87.00
A boa e velha ANOVA:
Com os dados acima, que mostram a densidade de uma praga em diferentes métodos de cultivo, nós vamos realizar uma análise de variância (carinhosamente conhecida por ANOVA) para testar a hipótese de que os métodos incluenciam a densidade da praga.
De novo, homogeneidade de variâncias
O mesmo pressuposto que vimos la no teste T se aplica aqui, e vamos testá-lo do mesmo jeito:
library(car) #Carregando o pacote necessário
leveneTest(pragas~cultivo)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.7813 0.5096
## 54Tudo beleza por aqui, pressuposto atendido!
E o pressuposto da normalidade?
Aqui nós vamos trabalhar com este pressuposto de uma maneira diferente. Veja só, a ANOVA faz parte de uma família de testes estatísticos chamados de modelos lineares. Nestes modelos, o fundamental é que os resíduos da análise sigam a distribuição normal. Então aqui as coisas mudam um pouco de ordem: primeiro criamos o modelo, depois verificamos o pressuposto.
Vamos lá, primeiro criando o modelo da ANOVA:
resultado<-aov(pragas~cultivo)Pois é, é só isso mesmo! Como salvamos o modelo dentro de um objeto, não podemos ver o resultado, a não ser que a gente o chame com a função summary. Mas não vamos fazer isso ainda; vamos, primeiro, conferir o tal pressuposto da normalidade dos resíduos, certo?
shapiro.test(resultado$residuals)##
## Shapiro-Wilk normality test
##
## data: resultado$residuals
## W = 0.9702, p-value = 0.1644Aqui tudo ok de novo! Então vamos conferir os resultados!
summary(resultado)## Df Sum Sq Mean Sq F value Pr(>F)
## cultivo 3 12322 4107 43.7 1.81e-14 ***
## Residuals 54 5076 94
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pronto!
A ANOVA nos indicou um resultado estatisticamente significativo. Neste caso, porém, isto quer dizer apenas que pelo menos duas médias diferem entre si. Se quisermos saber quais médias diferem de quem, temos que dar continuidade com a análise, e existem algumas abordagens distintas para se fazer isso. Vamos, por simplicidade, olhar para uma das mais comuns delas: os testes a posteirori.
Segue a execução do teste de Tukey:
TukeyHSD(resultado)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pragas ~ cultivo)
##
## $cultivo
## diff lwr upr p adj
## ausente-agroflorestal 41.1333333 31.582637 50.684030 0.0000000
## mata-agroflorestal 19.4666667 9.915970 29.017363 0.0000089
## refugio-agroflorestal 18.7857143 9.071746 28.499683 0.0000238
## mata-ausente -21.6666667 -31.051252 -12.282082 0.0000006
## refugio-ausente -22.3476190 -31.898316 -12.796922 0.0000005
## refugio-mata -0.6809524 -10.231649 8.869744 0.9975759Que nos indica que apenas um dos pares (cultivo próximo à mata e cultivo com áreas de refúgio) não diferem entre si.
Até aqui vimos duas funções novas:
aové a função que calcula a ANOVA;TukeyHSDrealiza o teste a posteriori de Tukey (e usa o resultado da anova como objeto);
Agora, que tal um gráfico?
Vamos fazer o mesmo tipo de gráfico que usamos para o teste T:
library(sciplot)
lineplot.CI(cultivo, pragas, type="p")
Bacana, mas notaram uma coisa: a classe “ausente”, que encaramos como um tipo de controle, ficou como a segunda classe na ordem do gráfico. Isto acontece porque o R usa ordem alfabética para organizar os níveis de um fator. Então precisamos mudar de ordem e refazer o gráfico. Podemos usar a função relevel para isso:
dados$cultivo<-relevel(cultivo, ref="ausente")
detach(dados)
attach(dados)
lineplot.CI(cultivo, pragas, type="p")
Agora sim!
Vamos mexer um pouco neste gráfico!
Vamos trabalhar mais neste último gráfico? Vamos imaginar que, por algum motivo, você queira que os nomes das categorias do eixo x sejam apresentadas na vertical, e não na horizontal (não há muito porque fazer isto neste exemplo, mas vamos seguir adiante apenas para explorar a questão, ok?).
Fazer isto é bem fácil, basta adicionar o argumento las=2, e pronto:
lineplot.CI(cultivo, pragas, type="p", las=2)
É, funcionou. Maaaaaaaas, ficou com probleminhas, não foi? Um dos nomes é grande para o espaço do gráfico, então temos que trabalhar nisso. Vamos, então, criar o gráfico com mais preparação e argumentos. E vamos aproveitar a deixa para mexer em mais coisas, ok? Veja como é, e vamos comentar tudo o que fizemos logo depois:
par(family="serif", mar=c(8, 4, 2, 2))
lineplot.CI(cultivo, pragas, type="p",
xlab="", ylab="Densidade da praga",
cex.lab=1.2, cex.axis=1.2, las=2, ylim=c(20, 80))
mtext("Tipo de cultivo", side=1, line=6.5,cex=1.2)
text(x=1, y=74, "a")
text(x=2, y=33, "c")
text(x=3, y=52, "b")
text(x=4, y=51, "b")
Ah, bem melhor! Vamos tentar entender tudo o que foi novidade aqui:
- o comando
parabre uma janela gráfica com os parâmetros que eu quiser estabelecer. Depois, todo o gráfico aberto atende àqueles parâmtros. - o argumento
familyserve para definir o tipo de fonte usada no gráfico; - o argumento
marserve para controlar as margens, os espaços em branco em volta da caixa do gráfico; - nós usamos
xlab=""para ele não escrever o nome do eixo, senão ele o faria em uma posição ruim, sobreposta aos nomes das classes; - a função
mtextadiciona texto na margem; neste caso, usamo-as para adicionar o nome do eixo x manualmente, na posição que escolhemo; - a função
textserve para se adicionar texto em qualquer canto do gráfico; no nosso caso, a usamos para adicionar letrinhas em cima da média, de acordo com o resultado do teste de Tukey;
Vamos fuçar um pouco na tal normalidade dos resíduos?
Apesar de termos feito o teste de normalidade, podemos desejar fazer uma inspeção visual dos nosso resíduos. Muitos autores, na verdade, defendem que o pressuposto da normalidade seja verificado assim, visualmente.
Usando um histograma
Um histograma pode ser usado para isso, apesar de não ser o melhor método para esta tarefa (veremos a maneira mais adequada a seguir, certo?). Mesmo assim, vamos ver como seria:
x<-scale(resultado$residuals)
hist(x, col="gray", main=NULL,
xlab="Resíduos", ylab="Frequências", las=1, freq=F,
ylim=c(0, 0.4))
curve(dnorm(x), col="red", lty=2, lwd=2, add=T)
Bonitinho, né? Vamos lá, detalhando algumas coisas que fizemos:
- a função
scalefoi usada para padronizar os dados dos resíduos, essencial para usarmos a curva sobreposta ao gráfico; - o argumento
freq=Ftambém é usado para permitir o uso da curva, e ele muda o eixo y para proporções, e não frequências; - a função
curveadiciona a curva ao histograma; - o argumento
add=Tgarante que a curva vai ser adicionada sem apagar o histogram; - o argumento
dnormcria o que seria esperado pela distribuição normal;
**Usando um “quantile-quantile plot**
O método mais adequado para se avaliar normalidade dos resíduos visualmente é criar o “quantile-quantile plot”, que veremos abaixo. A interpretação é bem tranquila: quanto mais os pontos estiverem seguindo a linha tracejada, mais eles se aproximam do que deveria ser esperado por uma distribuição normal.
Vejamos como fazer:
qqnorm(resultado$residuals)
qqline(resultado$residuals, lty=2)
Resumindo o que fizemos:
- a função
qqnormcria o plt, mas sem a linha tracejada; - a função
qqlineadiciona a linha tracejada no plot;
Ei, você disse que ia me ensinar o Kruskal-Wallis!
Ah é, já ia me esquecendo! :O
Em algumas situações, temos problemas de normalidade com os resíduos da ANOVA. Existem diversas saídas para se lidar com isso, e uma delas é o uso de um teste não paramétrico. O teste deste tipo que pode substitur a ANOVA é o Kruskal-Wallis.
No caso dos dados que estamos estudando, não há a menor necessidade de utilizá-lo. Mas vamos fazer isso assim mesmo só para ilustrar rapidamente os comandos necessários, ok?
O teste em si é super simples de ser executado:
kruskal.test(pragas~cultivo)##
## Kruskal-Wallis rank sum test
##
## data: pragas by cultivo
## Kruskal-Wallis chi-squared = 41.7865, df = 3, p-value = 4.453e-09E pronto, aí está o resultado.
Caso você esteja utilizando este teste justamente por conta de problemas de normalidade, pode ser legal representar seus resultados não com um gráfico de médias, e sim com um boxplot, que irá nos revelar um pouco mais sobre a estrutura dos dados:
boxplot(pragas~cultivo)
Só pela curiosidade, você pode comparar a visualização do boxplot com a do gráfico de médias. Aí a gente aproveita e aprende um argumentinho novo para ser usado dentro do par:
par(mfrow=c(1,2))
lineplot.CI(cultivo, pragas, type="p")
boxplot(pragas~cultivo)
- o argumento
mfrowdetermina que a janela de dados seja dividida em (linhas, colunas);
Ah, e se por acaso você precisar de um teste a posteriori para o Kruskal, há um disponível no pacote pgirmess. É um pacote grande e com várias dependências, então tem que ter paciência para baixar, ok? Aqui está o comando:
library(pgirmess)
kruskalmc(pragas~cultivo)## Multiple comparison test after Kruskal-Wallis
## p.value: 0.05
## Comparisons
## obs.dif critical.dif difference
## ausente-agroflorestal 40.4738095 16.55602 TRUE
## ausente-mata 20.9666667 16.26806 TRUE
## ausente-refugio 21.4380952 16.55602 TRUE
## agroflorestal-mata 19.5071429 16.55602 TRUE
## agroflorestal-refugio 19.0357143 16.83905 TRUE
## mata-refugio 0.4714286 16.55602 FALSE
E por hoje é só, espero que você tenha se divertido! ;)