O que é ANOVA e para que serve?

Os testes de hipóteses estão classificados em duas categorias: paramétricos e não paramétricos. Testes paramétricos são aqueles que utilizam os parâmetros da distribuição dos dados, ou uma estimativa destes, para o cálculo de sua estatística. Exemplo: Média, desvio-padrão e proporção. Já os testes não paramétricos não requerem tais pressupostos.

A Análise de Variância, ou ANOVA, é um teste paramétrico. Ela testa a hipótese de que a média de duas ou mais populações são iguais. A ANOVA é uma ferramenta que auxilia o pesquisador a avaliar a importância de um ou mais fatores, comparando as médias das variáveis resposta nos diferentes grupos.

Há duas formas de se utilizar a ANOVA:

O princípio utilizado pela ANOVA para determinar a diferença entre médias é baseada na análise de dois elementos da amostra: (i) a variação entre as médias dos grupos analisados; (ii) a variação em relação às amostras dentro do mesmo grupo.

Temos que:

SQ(total) = SQ(entre) + SQ(dentro)

Onde:

SQ(total) ou soma total de quadrados: é uma medida da variação total(em torno de x) em todos os dados amostrais combinados;

SQ(entre): é uma medida da variação entre as médias amostrais combinados.

Matematicamente calculamos o SQ(entre) da seguinte maneira:

Onde:

SQ(dentro) ou SQ(intra): é a soma de quadrados que representa a variabilidade comum a todas a populações em consideração.

Matematicamente calculamos o SQ(intra) da seguinte maneira:

Onde:

Após calculados as somas dos quadrados, calcularemos os graus de liberdade (GL) da nossa amostra da seguinte maneira:

GL(intra):

GL(entre):

Encontrados os nossos valores das somas dos quadrados, e os nossos graus de liberdade intra e entre grupos, podemos calcular os nossos quadrados médios (QM). Para isso, calculamos da seguinte maneira:

QM(intra) = SQ(intra)/GL(intra)
QM(entre) = SQ(entre)/GL(entre)

Após isso, podemos calcular o nosso valor de F, que diz respeito ao valor de poder que nos prediz quanto da fração é explicado pelo tratamento.

F é calculado da seguinte maneira:

Onde:

QM(entre): Variação entre as amostras

QM(intra): Variação dentro das amostras

É importante salientar que, para aplicarmos a ANOVA, os dados devem respeitar as seguintes premissas:

Lembrando, para o teste de hipótese usamos:

\(\mu\)A = \(\mu\)B = \(\mu\)C = \(\mu\)D

Caso o resultado seja estatisticamente significativo, aplica-se testes a posteriori para comparações múltiplas entre médias. Esses testes, chamados de post hoc, permitem identificar quais as populações são diferentes entre si, mantendo controlado o nível de significância do teste.

Delineamento amostral

Um pesquisador piauiense desenvolve junto ao Nuevo um estudo com bem-estar de macacos-prego (Sapajus libidinosus). Uma das suas perguntas de pesquisa é se há diferença nos níveis de estresse do macacos-prego nas diferentes áreas de seu trabalho, visto que o estresse afeta diretamente o bem-estar animal. Após consultar a literatura, ele decidiu utilizar os níveis de cortisol presente nas fezes dos macacos como variável preditora de estresse. Para isso, ele coletou as fezes de diferentes animais dos diferentes locais para verificar a quantidade de cortisol presente.

Configurando os dados para a ANOVA

Antes de importar os dados para o R, tenha certeza de que sua tabela está com a formatação correta. Aqui colocamos algumas dicas que vão te ajudar a deixar a sua planilha mais funcional:

Assim, uma planilha que continha alguns tipos de erros, como os destacados em vermelho da imagem a seguir:

Após corrigirmos os erros, a mesma deve ficar assim:

Após certificar-se que a planilha encontra-se apropriada, salve-a em formato .csv ou em formato .txt, para começarmos a trabalhar com os dados.

Usando seus dados no R

Após as duas etapas iniciais, importe seus dados para o R:

#Definindo diretório de trabalho

setwd("/home/daniel/Dados")

# Lendo a planilha

# Se estiver no formato .csv, use

macacos <- read.csv("dados_macacos.csv")

# Se estiver no formato .txt, use

macacos2 <- read.delim("dados_macacos.txt")

# OBS: lembre-se de especificar o argumento sep = ";"

# Caso a configuração do seu computador esteja com separador decimal como ",", chame a tabela da seguinte maneira no R:

macacos3 <- read.csv2("dados_macacos.csv", dec = ",")

# Visualize sua planilha de dados

View(macacos)

Aplicando a ANOVA

A função do R que executa a ANOVA é a aov. Nessa função temos como argumentos nossa a variável resposta(cortisol), variável preditora (local) e nossa planilha de dados:

model <- aov(Cortisol ~ Local, data = macacos)
model
## Call:
##    aov(formula = Cortisol ~ Local, data = macacos)
## 
## Terms:
##                   Local Residuals
## Sum of Squares  6248447    324986
## Deg. of Freedom       3       116
## 
## Residual standard error: 52.93016
## Estimated effects may be unbalanced

Essa função nos retorna alguns elementos:

Com a função summary conseguimos mais informações:

summary(model)
##              Df  Sum Sq Mean Sq F value Pr(>F)    
## Local         3 6248447 2082816   743.4 <2e-16 ***
## Residuals   116  324986    2802                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Temos agora informação de:

A partir da estatística F e seu valor-p abaixo de 0.05 podemos temos embasamento estatístico para afirmar com grande confiança que as médias de Cortisol dos macacos difere entre os locais analisado.

Testando as premissas da ANOVA

Embora tenhamos encontrado que há um diferença entre as médias, o resultado da ANOVA só é robusto se as premissas do testes forem satisfeitas.

Homogeneidade das amostras

Uma vez que a ANOVA tem como premissa a homogeneidade dos dados, podemos testar se nossos dados são adequados nesse quesito com o Teste de Levene para homocedasticia. Para isso precisaremos do pacote car:

library(car)
leveneTest(Cortisol ~ Local, data = macacos)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   3  0.6416 0.5898
##       116

A hipótese nula do Teste de Levene é de que não há diferença entre as variâncias dos grupos. O valor-p maior do que 0.05 nos dá uma confiança estatística para afirmar que as variâncias são de fato iguais e portanto nossos dados são homogêneos.

Normalidade dos resíduos

Já premissa da ANOVA referente a normalidade dos resíduos pode ser testada através do teste de Shapiro-Wilk:

shapiro.test(resid(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98995, p-value = 0.5294

A hipótese nula do Teste de Shapiro-Wilk é de que não há diferença entre a nossa distribuição dos dados e a distribuição normal. O valor-p maior do que 0.05 nos dá uma confiança estatística para afirmar que as distribuição dos nossos resíduos não difere da distribuição normal.

Dessa forma nossos dados satisfazem todas as premissas da ANOVA e portanto, o resultado da nossa ANOVA são válidos.

Comparação de médias par-a-par

Quando fizemos a ANOVA, rejeitamos a hipótese nula de que todas as médias são iguais, no entanto, não sabemos exatamente quais grupos têm médias diferentes. Para avaliar isso podemos fazer um Teste de Tukey:

TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cortisol ~ Local, data = macacos)
## 
## $Local
##                                 diff        lwr        upr     p adj
## Pedra Furada-Baixa Grande -154.16410 -189.78811 -118.54009 0.0000000
## Recife-Baixa Grande        457.59546  421.97145  493.21948 0.0000000
## Teresina-Baixa Grande       13.11672  -22.50729   48.74073 0.7724681
## Recife-Pedra Furada        611.75957  576.13556  647.38358 0.0000000
## Teresina-Pedra Furada      167.28083  131.65682  202.90484 0.0000000
## Teresina-Recife           -444.47874 -480.10275 -408.85473 0.0000000

A função retorna uma tabela onde as linhas são referentes à comparação par-a-par entre os grupos:

No nosso exemplo temos que todas as médias de cortisol diferem entre todos os locais, exceto entre Teresina e Baixa Grande.

Gráficos

É importante também exibirmos graficamente a distribuição dos nossos dados, uma vez que esta é uma ferramenta que facilita o entendimento do comportamento dos dados.

Para produzir os gráficos a seguir utilizamos os pacotes ggplot2 e ggpubr.

Histograma

Os gráficos de histograma nos mostram como os dados estão distribuídos. Com esse tipo de gráfico é possível estimar onde os valores estão concentrados, quais são os extremos e se existem lacunas ou valores incomuns. Eles também são úteis ao fornecer uma visão aproximada da distribuição de probabilidade.

library(ggplot2)
ggplot(macacos, aes(x = Cortisol,fill = Local)) +
  geom_histogram(color = "black", binwidth = 50)+
  facet_grid(Local ~ .) +
  labs(y = 'Frequência') +
  scale_fill_manual(values=c("#0f8bf7", "#f7830f", "#ff12d0", "#2a8008"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = 'top',
        axis.line = element_line(colour = "black"))

O eixo x indica o nível de cortisol. O eixo y indica a frequência das observações.

Boxplot

O boxplot, assim como o histograma, mostra a distribuição do dados, porém exibindo-a através de quartis. Os 3 segmentos que formam o “box” equivalem ao 1º quartil, 2º quartil (ou mediana) e 3º quartil. As linhas verticais mostram os limites mínimos e máximos dos valores da amostra. O boxplot também exibe os outliers (representados pelos pontos), que são amostras que divergem muito da distribuição do restante das amostras.

library(ggpubr)
ggboxplot(macacos, x = "Local", y = "Cortisol",
          fill = "Local", palette = c("#0f8bf7", "#f7830f", "#ff12d0", "#2a8008"),
          order = c("Baixa Grande","Pedra Furada", "Recife", "Teresina"),
          ylab = "Níveis de cortisol", xlab = "Locais de estudo")

O eixo x indica os locais. O eixo y indica o nível de cortisol.

Gráfico Violino

O gráfico de violino, ou como também é chamado beanplot, é similar ao boxplot. É usado para visualizar a distribuição dos dados e a densidade de probabilidade. O gráfico é uma combinação do boxplot e gráfico de densidade, que é girado na vertical e espelhado, para mostrar a distribuição dos dados.

library(ggpubr)
ggviolin (macacos, x = "Local", y = "Cortisol", fill = "Local",
          palette = c("#0f8bf7", "#f7830f", "#ff12d0", "#2a8008"),
          order = c("Baixa Grande","Pedra Furada", "Recife", "Teresina"),
          add = "boxplot", add.params = list(fill = "white"),
          ylab = "Níveis de cortisol", xlab = "Locais de estudo")

O eixo x indica os locais. O eixo y indica o nível de cortisol.

Barplot

O barplot representa as médias dos grupos e seu respectivo intervalo de confiança

#Calculando o intervalo de confiança da média
fatores <- unique(macacos$Local)
data.bp <- data.frame(matrix(nrow = length(fatores), ncol = 4))
for (i in 1:length(fatores)){
  media <- mean(macacos$Cortisol[macacos$Local == fatores[i]])
  sd <- sd(macacos$Cortisol[macacos$Local == fatores[i]])
  n <- length(macacos$Cortisol[macacos$Local == fatores[i]])
  se <- sd(macacos$Cortisol[macacos$Local == fatores[i]]) /
  sqrt(length(macacos$Cortisol[macacos$Local == fatores[i]]))
  alpha = 0.05
  t = qt((1 - alpha) / 2 + 0.5, length(macacos$Cortisol[macacos$Local == fatores[i]]) - 1)
  ci = t * se
  lower <- media - ci
  upper <- media + ci
  data.bp[i, 2:4] <- data.frame(media, lower, upper)
}

data.bp[, 1] <- fatores
names(data.bp) <- c('Local', 'media', 'lower', 'upper')
head(data.bp)
##          Local    media    lower    upper
## 1 Baixa Grande 452.4705 432.9370 472.0040
## 2 Pedra Furada 298.3064 277.9987 318.6141
## 3     Teresina 465.5872 448.2114 482.9630
## 4       Recife 910.0660 888.4637 931.6682
#Plotando o barplot
library(ggplot2)
ggplot(data.bp) +
  geom_bar(aes(x = Local, y = media, fill = Local), stat = "identity") +
 
  geom_errorbar( aes(x = Local, ymin = lower, ymax= upper),
                 width = 0.5, colour = "black") +
  scale_fill_manual(values=c("#0f8bf7", "#f7830f", "#ff12d0", "#2a8008")) +
  labs(x = 'Local', y = 'Média de Cortisol') +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = 'top',
        axis.line = element_line(colour = "black"))

O eixo x indica os locais. O eixo y as médias de cortisol na fezes dos macacos em cada um dos locais. As linha horizontais acima das barras indicam os valores mínimos e máximos do intervalo de confiança da média de cortisol.