ANOVA e teste de Tukey
Análise de variância (ANOVA) e testes de médias são métodos comuns em artigos científicos. Você com certeza já viu aquelas letrinhas indicando a diferença entre tratamentos em algum estudo publicado. Por mais que este método esteja entrando em desuso - há uma tendência em abandonar esse tipo de abordagem estatística - penso que ainda o veremos por muitos anos no meio científico.
Como contexto, temos um teste de 5 progênies de eucalipto e queremos avaliar se volume por hectare (nossa variável resposta), difere entre os tratamentos.
Pois bem, para percebermos a dimensão dos dados e qual a variabilidade de cada tratamento, vamos criar um boxplot (Figura 1). Caso você queira saber um pouco mais sobre este tipo de gráfico, veja o post sobre ele.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readr, dplyr, tibble, ggplot2, car, agricolae)dados <- read_csv2(
"https://github.com/italocegatta/italocegatta.github.io_source/raw/master/content/dados/base_progenie.csv"
)
dados## # A tibble: 30 x 3
## repeticao progenie volume
## <dbl> <chr> <dbl>
## 1 1 A 212
## 2 2 A 206
## 3 3 A 224
## 4 4 A 289
## 5 5 A 324
## 6 6 A 219
## 7 1 B 108
## 8 2 B 194
## 9 3 B 163
## 10 4 B 111
## # ... with 20 more rowsggplot(dados, aes(progenie, volume)) +
geom_boxplot(fill = "grey60", alpha = 0.8) +
theme_bw(16)
Figura 1: Variabilidade do volume por hectare de cada tratamento.
A ANOVA é um método bastante consolidado no meio acadêmico. Basicamente, este método informa se existe um tratamento discrepante dentre os demais. Entretanto, ele exige que algumas premissas sejam atendidas, como: distribuição normal dos resíduos e homogeneidade de variância.
Primeiro, vamos utilizar o teste de Levene para verificar se há homogeneidade de variância, ou homocedasticidade. Como o p-valor é maior que 5% não temos evidência significativa para rejeitar a hipótese nula de homogeneidade, ou seja, nossos dados tem homogeneidade de variância.
leveneTest(volume ~ factor(progenie), data=dados)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 2.4677 0.07086 .
## 25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O segundo pressuposto é a normalidade dos resíduos. Utilizaremos o teste de Shapiro-Wilk cuja hipótese nula é a de que os dados seguem uma distribuição normal. Como o p-valor é superior ao limite de 5%, podemos aceitar a hipótese nula e considerar nossos dados normais.
anova <- aov(volume ~ progenie, data=dados)
shapiro.test(resid(anova))##
## Shapiro-Wilk normality test
##
## data: resid(anova)
## W = 0.96097, p-value = 0.3279Uma vez que os pressupostos foram atendidos, seguiremos para a ANOVA. Note que, caso os testes de Levene e Shapiro-Wilk resultassem em um p-valor significante, ou seja, menor que 5%, teríamos que utilizar outro método estatístico para analisar nossos dados. Nesse caso, uma alternativa é utilizar testes não-paramétricos, uma vez que eles não exigem os pressupostos que acabamos de testar.
Nossa ANOVA resultou em um p-valor menor que 5%, portanto, temos evidências de que ao menos um tratamento se diferencia dos demais. Isso já é uma resposta, mas pouco acrescenta à nossa pesquisa pois queremos saber quem é este tratamento discrepante. Ou melhor, queremos poder comparar os tratamentos entre si e verificar quais são estatisticamente iguais ou diferentes.
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## progenie 4 86726 21681 8.89 0.000131 ***
## Residuals 25 60974 2439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para esta abordagem existem alguns testes de médias e cada um tem uma particularidade, mas de longe o mais utilizado é o de Tukey.
A interpretação do teste de Tukey é simples. Após determinarmos a diferença mínima significativa (ou Honest Significant Difference - HSD), podemos julgar se as médias são iguais ou não. Em termos práticos, esse valor nos dá uma margem de igualdade, pois se a diferença entre dois tratamentos for maior do que isso, os médias são diferentes.
A análise começa sempre pela maior média, no nosso caso a progênie A (245, 66). Com uma continha rápida, a média do tratamento A menos a diferença mínima significativa 245,66 - 83,73 = 161,93, aceitaremos que um tratamento é igual ao A se a média dele for maior que 161,93. O tratamento subsequente (o segundo do ranking) é a progênie D e como sua média é maior que 161,93 podemos dizer que ela é estatisticamente igual a progênie A.
As próximas comparações seguem a mesma lógica. Quando registramos que duas médias são iguais, nós as rotulamos com a mesma letra para facilitar a identificação. Veja no fim do output as letras evidenciando a igualdade entre os tratamentos.
tukey <- HSD.test(anova, "progenie")
tukey## $statistics
## MSerror Df Mean CV MSD
## 2438.953 25 165.7667 29.79233 83.73866
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey progenie 5 4.153363 0.05
##
## $means
## volume std r Min Max Q25 Q50 Q75
## A 245.6667 48.78798 6 206 324 213.75 221.5 272.75
## B 159.6667 49.47996 6 108 236 119.75 154.5 186.25
## C 80.5000 15.60449 6 63 100 70.00 76.5 93.50
## D 190.1667 75.37484 6 100 267 121.75 207.0 251.75
## E 152.8333 37.96534 6 106 210 133.75 141.5 175.50
##
## $comparison
## NULL
##
## $groups
## volume groups
## A 245.6667 a
## D 190.1667 ab
## B 159.6667 bc
## E 152.8333 bc
## C 80.5000 c
##
## attr(,"class")
## [1] "group"Para deixar mais visual ainda, podemos construir um gráfico de barras com a média de cada tratamento e adicionar a sua letra correspondente ao teste de Tukey (Figura 2).
tukey$groups %>%
rownames_to_column(var = "trt") %>%
mutate(trt = reorder(trt, -volume, mean)) %>%
ggplot(aes(trt, volume)) +
geom_col(alpha = 0.8, color = "black") +
geom_text(aes(label = groups), vjust = 1.8, size = 9, color = "white") +
labs(x = "Progênies", y = "Médias") +
theme_bw(16)
Figura 2: Médias dos tratamentos. As letras indicam médias estatisticamente iguais pelo teste de Tukey a 5% de significância.
Caso tenha alguma dúvida ou sugestão sobre o post, fique à vontade para fazer um comentário ou me contatar por E-mail.
sessioninfo::session_info(c("readr", "dplyr", "tibble", "ggplot2", "car", "agricolae"))## - Session info ----------------------------------------------------------
## setting value
## version R version 3.5.3 (2019-03-11)
## os Windows 10 x64
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate Portuguese_Brazil.1252
## ctype Portuguese_Brazil.1252
## tz America/Sao_Paulo
## date 2019-08-25
##
## - Packages --------------------------------------------------------------
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 3.5.0)
## agricolae * 1.3-1 2019-04-04 [1] CRAN (R 3.5.3)
## AlgDesign 1.1-7.3 2014-10-15 [1] CRAN (R 3.5.2)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.5.3)
## backports 1.1.4 2019-04-10 [1] CRAN (R 3.5.3)
## BH 1.69.0-1 2019-01-07 [1] CRAN (R 3.5.2)
## boot 1.3-20 2017-08-06 [2] CRAN (R 3.5.3)
## car * 3.0-3 2019-05-27 [1] CRAN (R 3.5.3)
## carData * 3.0-2 2018-09-30 [1] CRAN (R 3.5.2)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.5.1)
## class 7.3-15 2019-01-01 [2] CRAN (R 3.5.3)
## classInt 0.4-1 2019-08-06 [1] CRAN (R 3.5.3)
## cli 1.1.0 2019-03-19 [1] CRAN (R 3.5.3)
## clipr 0.7.0 2019-07-23 [1] CRAN (R 3.5.3)
## cluster 2.0.7-1 2018-04-13 [2] CRAN (R 3.5.3)
## coda 0.19-3 2019-07-05 [1] CRAN (R 3.5.3)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.5.3)
## combinat 0.0-8 2012-10-29 [1] CRAN (R 3.5.2)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
## curl 4.0 2019-07-22 [1] CRAN (R 3.5.3)
## data.table 1.12.2 2019-04-07 [1] CRAN (R 3.5.3)
## DBI 1.0.0 2018-05-02 [1] CRAN (R 3.5.1)
## deldir 0.1-23 2019-07-31 [1] CRAN (R 3.5.3)
## digest 0.6.20 2019-07-04 [1] CRAN (R 3.5.3)
## dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.5.3)
## e1071 1.7-2 2019-06-05 [1] CRAN (R 3.5.3)
## ellipsis 0.2.0.1 2019-07-02 [1] CRAN (R 3.5.3)
## expm 0.999-4 2019-03-21 [1] CRAN (R 3.5.3)
## fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.1)
## forcats 0.4.0 2019-02-17 [1] CRAN (R 3.5.2)
## foreign 0.8-71 2018-07-20 [2] CRAN (R 3.5.3)
## gdata 2.18.0 2017-06-06 [1] CRAN (R 3.5.2)
## ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.5.3)
## glue 1.3.1 2019-03-12 [1] CRAN (R 3.5.3)
## gmodels 2.18.1 2018-06-25 [1] CRAN (R 3.5.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.5.3)
## gtools 3.8.1 2018-06-26 [1] CRAN (R 3.5.2)
## haven 2.1.1 2019-07-04 [1] CRAN (R 3.5.3)
## highr 0.8 2019-03-20 [1] CRAN (R 3.5.3)
## hms 0.5.0 2019-07-09 [1] CRAN (R 3.5.3)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.1)
## httpuv 1.5.1 2019-04-05 [1] CRAN (R 3.5.3)
## jsonlite 1.6 2018-12-07 [1] CRAN (R 3.5.2)
## KernSmooth 2.23-15 2015-06-29 [2] CRAN (R 3.5.3)
## klaR 0.6-14 2018-03-19 [1] CRAN (R 3.5.2)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
## labelled 2.2.1 2019-05-26 [1] CRAN (R 3.5.3)
## later 0.8.0 2019-02-11 [1] CRAN (R 3.5.2)
## lattice 0.20-38 2018-11-04 [2] CRAN (R 3.5.3)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.5.3)
## LearnBayes 2.15.1 2018-03-18 [1] CRAN (R 3.5.2)
## lme4 1.1-21 2019-03-05 [1] CRAN (R 3.5.3)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
## maptools 0.9-5 2019-02-18 [1] CRAN (R 3.5.2)
## MASS 7.3-51.1 2018-11-01 [2] CRAN (R 3.5.3)
## Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.5.3)
## MatrixModels 0.4-1 2015-08-22 [1] CRAN (R 3.5.2)
## mgcv 1.8-28 2019-03-21 [1] CRAN (R 3.5.3)
## mime 0.7 2019-06-11 [1] CRAN (R 3.5.3)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 3.5.2)
## minqa 1.2.4 2014-10-09 [1] CRAN (R 3.5.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
## nlme 3.1-137 2018-04-07 [2] CRAN (R 3.5.3)
## nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.5.1)
## nnet 7.3-12 2016-02-02 [2] CRAN (R 3.5.3)
## openxlsx 4.1.0.1 2019-05-28 [1] CRAN (R 3.5.3)
## pbkrtest 0.4-7 2017-03-15 [1] CRAN (R 3.5.2)
## pillar 1.4.2 2019-06-29 [1] CRAN (R 3.5.3)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
## plogr 0.2.0 2018-03-25 [1] CRAN (R 3.5.1)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1)
## progress 1.2.2 2019-05-16 [1] CRAN (R 3.5.3)
## promises 1.0.1 2018-04-13 [1] CRAN (R 3.5.1)
## purrr 0.3.2 2019-03-15 [1] CRAN (R 3.5.3)
## quantreg 5.51 2019-08-07 [1] CRAN (R 3.5.3)
## questionr 0.7.0 2018-11-26 [1] CRAN (R 3.5.2)
## R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.2)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.5.0)
## Rcpp 1.0.2 2019-07-25 [1] CRAN (R 3.5.3)
## RcppEigen 0.3.3.5.0 2018-11-24 [1] CRAN (R 3.5.1)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.5.2)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 3.5.3)
## rematch 1.0.1 2016-04-21 [1] CRAN (R 3.5.1)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.1)
## rio 0.5.16 2018-11-26 [1] CRAN (R 3.5.2)
## rlang 0.4.0 2019-06-25 [1] CRAN (R 3.5.3)
## rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.5.3)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
## sf 0.7-7 2019-07-24 [1] CRAN (R 3.5.3)
## shiny 1.3.2 2019-04-22 [1] CRAN (R 3.5.3)
## sourcetools 0.1.7 2018-04-25 [1] CRAN (R 3.5.1)
## sp 1.3-1 2018-06-05 [1] CRAN (R 3.5.1)
## SparseM 1.77 2017-04-23 [1] CRAN (R 3.5.2)
## spData 0.3.0 2019-01-07 [1] CRAN (R 3.5.2)
## spdep 1.1-2 2019-04-05 [1] CRAN (R 3.5.3)
## stringi 1.4.3 2019-03-12 [1] CRAN (R 3.5.3)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 3.5.2)
## tibble * 2.1.3 2019-06-06 [1] CRAN (R 3.5.3)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
## units 0.6-3 2019-05-03 [1] CRAN (R 3.5.3)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 3.5.1)
## vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.5.3)
## viridisLite 0.3.0 2018-02-01 [1] CRAN (R 3.5.1)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 3.5.3)
## zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.5.2)
## zip 2.0.3 2019-07-03 [1] CRAN (R 3.5.3)
##
## [1] C:/Users/Italo/Documents/R/win-library/3.5
## [2] C:/Program Files/R/R-3.5.3/library