Estatística/Análise de Dados


Análise de Variância (ANOVA)


Prof. Rodrigo Sant'Ana

UNIVALI / CTTMar / GEP


Contextualização


A Análise de Variância (ANOVA) trata-se de um método estatı́stico que permite realizar comparações simultâneas entre duas ou mais médias, ou seja, permite testar hipóteses sobre médias de distintas populações.

Em um caso simples, onde temos uma variável categórica X (e.g. variável explicativa do modelo) que contenha mais de duas categorias (exemplo: estações do ano, k = 4 categorias - verão, outono, inverno e primavera), as medições de uma dada variável resposta Y de nosso interesse passa a ser induzida à um particionamento em subpopulações, neste caso 4 subpopulações - estações do ano. E é sobre estas subpopulações que construímos a hipótese a ser testada.

As diferenças eventuais entre cada das Xi subpopulações (com i = k número de estações do ano) poderiam ser testadas tomando-as duas a duas e utilizando testes de diferenças entre dois grupos (exemplo: teste t). No entanto, a combinação destes procedimentos dois à dois acarretaria em um número elevado de comparações, onde este quantitativo pode ser estimado a priori pela equação:


$$C_{k}^{2}=\frac{k*(k-1)}{2}$$

Mesmo que o problema não esteja fadado ao número de testes a serem realizados, há uma implicação muito mais complexa que este esforço em realizar testes combinados, relacionada diretamente ao controle do nível de significância desta combinação. De maneira geral, a probabilidade de erros do tipo I e tipo II é acumulada a cada combinação, de maneira que, quando executamos este tipo de combinação a probabilidade de termos um erro do tipo I pode ser estimada por:


αcombinação = 1 − (1 − α)c


onde α é o nível de significância para cada teste t realizado e c é o número de testes t independentes. Ou seja, em um caso onde temos quatro subpopulações a serem testadas, teriamos que realizar 6 testes t independentes, resultantes da combinação entre as 4 estações do ano (verão x inverno, verão x primavera, verão x outono, inverno x primavera, inverno x outono, primavera x outono). Considerando um nível de significância para cada teste independente α = 0.05, a probabilidade real de termos um erro do Tipo I não seria igual a 5% como esperado inicialmente, mas sim de 26,5% (0,265) estimado pela combinação dos níveis de significância dos testes independentes.

Neste sentido, a Análise de Variância (ANOVA) se justifica, pois permite testar a diferença entre os k grupos globalmente, utilizando um único teste estatístico.


ANOVA


Análise de Variância trata-se de um teste sobre a igualdade (ou não) dos valores esperados (médias) de uma determinada variável de interesse nas k subpopulações de interesse. Assim, o que se desejar testar é:


H0 : μ1 = μ2 = ... = μk
H1 :  pelo menos um μ é distinto dos demais.

O teste é construído comparando-se duas estimativas independentes de σ2 a partir das k subpopulações. A primeira expressa eventuais diferenças entre as subpopulações (ou grupos), mais precisamente sob os desvios entre as estimativas médias das subpopulações (μ1, μ2, ..., μk) e a estimativa da média populacional (μ) (SQEnt - Soma dos Quadrados entre os Grupos). Já a segunda, esta centrada na diferença dentro dos grupos, focada nos desvios entre as observações e a média amostral de seu respectivo grupo (SQDen - Soma dos Quadrados dentro dos Grupos).

Assim, considerando este particionamento da variabilidade total (SQT - Soma dos Quadrados Totais) e seus respectivos graus de liberdade, pode-se calcular duas estimativas não viciadas da variância populacional σ2. Uma baseada na razão entre SQDen e seus respectivos graus de liberdade (n − k) e, outra, dada pela razão entre SQEnt e seus graus de liberdade (k − 1).

Por fim, o teste se motiva em comparar estas duas medidas da variância populacional não viciada, onde a razão entre elas segue uma distribuição F com k − 1 e n − k graus de liberdade. Se pensarmos algebricamente, o que esperaríamos desta razão, caso H0 fosse verdadeira, é que o valor resultante deste quociente fosse próximo ou igual a 1, assumindo assim que ambas estimações da variância populacional são iguais.

Pressupostos ou Suposições do Método

Assim como outros testes de hipóteses, a Análise de Variância também se estrutura sob algumas suposições ou pressupostos para que seja aplicável, sendo estes:

  1. Todas as observações devem ser independentes;
  2. As observações em cada grupo devem possuir uma distribuição, aproximadamente normal;
  3. As variâncias em cada grupo devem ser aproximadamente iguais.


A independência entre as observações é sempre importante em uma ANOVA. A condição de normalidade é muito importante quando se têm pequenas amostras em cada grupo. Já a condição de constância das variâncias é especialmente importante quando os tamanhos das amostras que se pretende analisar são diferentes.

(Çetinkaya-Rundel, 2015 - Statistical Analysis and Inference Course / Diez et al. 2015, OpenIntro Statistics)



ANOVA no R

Análise de Variância Simples

A base de dados abaixo consiste em amostras de intensidade de vento provenientes de 4 áreas distintas onde serão instalados parques eólicos. No entanto, a empresa responsável pela instalação está entrando no mercado nacional agora, e pouco conhece a realidade destas áreas, tampouco a constância de vento em cada uma delas. E com isso, está um pouco receosa se vale a pena concorrer nos 4 processos licitatórios de cada uma das áreas.

Para que a empresa não tome uma decisão às escuras, foram disponibilizados dados de medições de vento em cada uma das regiões. A ideia é avaliar se nas quatro regiões as intensidades de vento são iguais. Caso contrário, identificar qual ou quais regiões não vale a pena participar do processo de licitação.

As informações disponibilizadas seguem apresentadas na tabela abaixo. Todos os dados foram disponibilizados em metros por segundo (m/s).

Tabela 1 - Medições de intensidade de vento nas 4 (quatro) áreas onde serão instalados os parques eólicos.

ID Área I Área II Área III Área IV
1 3,2 4,2 5,4 4,5
2 3,5 3,7 4,6 3,8
3 2,7 3,4 4 4,1
4 4,1 4,3 5,3 3,1
5 3,1 3,9 4,7 4,2
6 3,7 4,1 4,2 3,4
7 4,2 3,1 4,9 4,2
8 3,6 4,5 4,7 4,5

Com o intuito de responder está questão, avalie se os pressupostos da Análise de Variância são atendidos, construa a hipóstese a ser testada à um nível de significância α = 0, 05, aplique a ANOVA e interprete os resultados alcançados.

Preparando o RStudio para o Trabalho

  1. Criar uma pasta de trabalho no seu computador onde serão armazenadas todas as informações, análises e documentos que deverão ser utilizados durante esta pesquisa/trabalho;

  2. Digite os dados em uma planilha eletrônica, respeitando os conceitos de organização de dados;

  3. Salve esta planilha no formato csv, dentro da pasta de trabalho criada anteriormente;

  4. Abra o RStudio;

  5. Direcione o RStudio para sua pasta de trabalho criada para esta aula;

  6. No RStudio, crie um novo script para receber o algoritmo (comandos) que sustentarão as análises deste trabalho;

  7. Salve o script na pasta de trabalho criada anteriormente;

  8. Importe a base de dados para o R utilizando os conhecimentos aprendidos em sala;

  9. Feitos os passos acima, é hora de começar as análises.

Realizando a análise no R

  • Importando a base de dados...
R source
### Carregando os dados no R utilizando a importação a partir de um
### arquivo salvo no computador
dados <- read.table("vento.csv", header = TRUE, sep = ";", dec = ",")
  • Verificando o pressuposto de homogeneidade das variâncias...

Uma maneira simples é visualizar as variâncias para cada um dos tratamentos e verificar se elas distoam muito entre si. Neste exemplo podemos notar que as variâncias amostrais calculadas para cada uma das regiões são bastante próximas.

R source
### Visualizando as variâncias de cada uma das regiões
aggregate(vento ~ area, dados, var)
R output
  area     vento
1    I 0.2555357
2   II 0.2257143
3  III 0.2335714
4   IV 0.2564286

Outra maneira é testar a hipótese de homogeneidade das variâncias através de um teste específico - exemplo o teste de Bartlett. Notem o p-valor do teste, ele resume o resultado de um teste de hipótese. Quando p-valor for menor que o nível de significância escolhido para o teste (em geral, convencionado em 0,05), rejeita-se a hipótese nula em favorecimento da hipótese alternativa - neste caso assume-se a desigualdade das variâncias; Caso ele seja maior que o nível de significância, falha-se em rejeitar H0 por falta de evidências, sendo assim, conclui-se que as variâncias são constantes ou homogêneas.

No caso do exemplo dos ventos, o p-valor do teste de Bartlett foi maior que 0,05, assim há homogeneidade das variâncias.

R source
### Testando a homogeneidade das variâncias através do teste de
### Bartlett. 
bartlett.test(vento ~ area, dados)
R output

    Bartlett test of homogeneity of variances

data:  vento by area
Bartlett's K-squared = 0.040882, df = 3, p-value = 0.9978
  • Verificando a normalidade dos dados...

Uma alternativa para verificar o pressuposto de normalidade é o teste de Shapiro-Wilk. Notem que a percepção sobre o resultado do teste também se estrutura em torno do p-valor do mesmo modo que visto para o teste de homogeneidade das variâncias.

H0: consistem em dizer que há normalidade nos dados; H1: que não há.

Assim, rejeita-se H0 quando p − valor < α, por convenção utiliza-se α =  0,05.

Neste exemplo, os dados de vento também atendem o pressuposto de normalidade da ANOVA (p-valor > 0,05).

R source
### Testando a normalidade dos dados através do teste de Shapiro-Wilk
shapiro.test(dados$vento[dados$area == "I"])
R output

    Shapiro-Wilk normality test

data:  dados$vento[dados$area == "I"]
W = 0.96981, p-value = 0.8965
R source
shapiro.test(dados$vento[dados$area == "II"])
R output

    Shapiro-Wilk normality test

data:  dados$vento[dados$area == "II"]
W = 0.95848, p-value = 0.7955
R source
shapiro.test(dados$vento[dados$area == "III"])
R output

    Shapiro-Wilk normality test

data:  dados$vento[dados$area == "III"]
W = 0.95192, p-value = 0.7305
R source
shapiro.test(dados$vento[dados$area == "IV"])
R output

    Shapiro-Wilk normality test

data:  dados$vento[dados$area == "IV"]
W = 0.89956, p-value = 0.2863
  • Aplicando a ANOVA...

No R, existem duas funções básicas capazes de aplicar a Análise de Variância, a função anova() e a função aov().

R source
### Análise de Variância no R - Forma I
tab.anova <- anova(lm(dados$vento ~ dados$area))
tab.anova
R output
Analysis of Variance Table

Response: dados$vento
           Df Sum Sq Mean Sq F value    Pr(>F)    
dados$area  3 6.1659 2.05531  8.4646 0.0003682 ***
Residuals  28 6.7987 0.24281                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R source
### Análise de Variância no R - Forma II
tab.anova2 <- aov(dados$vento ~ dados$area)
summary(tab.anova2)
R output
            Df Sum Sq Mean Sq F value   Pr(>F)    
dados$area   3  6.166  2.0553   8.465 0.000368 ***
Residuals   28  6.799  0.2428                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Em ambas saídas da análise, pode-se notar a existência de uma diferença significativa entre os padrões de ventos das quatro regiões. Notem que o p-valor estimado foi igual à 0,000368 (Pr(>F)). Isto nos leva a rejeitar a hipótese nula em favorecimento à hipótese alternativa, ou seja, pelo menos uma das regiões apresenta um padrão médio de vento distinto das demais.

Antes de verificar qual ou quais regiões são distintas entre sí, é importante verificar o pressuposto de independência das observações. Este pressuposto pode ser avaliado através da análise dos resíduos do modelo ANOVA.

Aplicando a função genérica plot() do R aos resíduos do modelo, que podem ser acessados pela função residuals() ao objeto que armazena os resultados do modelo ANOVA pode-se visualizar a distribuição dos resíduos e sua independência à cada observação. Outra forma de visualizar esta independência dos dados sobre a análise dos resíduos de um modelo esta centrada na distribuição dos mesmos. Espera-se que os resíduos de um modelo bem ajustado possua uma distribuição normal, centrada em 0 e variância constante.

R source
### Para verificar a dependência nas observações, podemos analisar os
### residuos do modelo ANOVA. 
par(mfrow = c(1, 2))
plot(residuals(tab.anova2), xlab = "Indice da observação",
     ylab = "Resíduos", ylim = c(-1, 1))
hist(residuals(tab.anova2), include.lowest = TRUE, right = FALSE,
     main = "", xlab = "Resíduos", ylab = "Frequência",
     col = "gray")
R source
par(mfrow = c(1, 1))

Outro resultado importante de um modelo pode ser acessado pelo cálculo do coeficiente de explicação r2. O coeficiente de explicação de um modelo ANOVA é dados pela razão entre a Soma dos Quadrados dos Resíduos entre os Grupos (SQEnt)e a Soma dos Quadrados dos Resíduos Totais (SQT). Utilizando a saída da função_anova()_ fica mais simples acessar os valores para estimar o r2.

R source
### calculando o coeficiente de explicação - r2
tab.anova$"Sum Sq"[1]/sum(tab.anova$"Sum Sq")
R output
[1] 0.4755948

Por fim, resta identificar qual ou quais áreas possuem média de ventos distintas entre sí. Para isto, utilizamos a função TukeyHSD() do R para contrastar as diferenças. Note que esta função deve ser aplicada sobre o modelo ANOVA ajustado com a função aov().

R source
### Verificando quais médias são significativamente distintas entre si,
### através da aplicação de um teste a Post-Hoc ou Posterior, chamado
### Teste de Tukey
post.hoc <- TukeyHSD(tab.anova2)
post.hoc
R output
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = dados$vento ~ dados$area)

$`dados$area`
          diff        lwr         upr     p adj
II-I    0.3875 -0.2851943  1.06019428 0.4098324
III-I   1.2125  0.5398057  1.88519428 0.0001921
IV-I    0.4625 -0.2101943  1.13519428 0.2604982
III-II  0.8250  0.1523057  1.49769428 0.0117845
IV-II   0.0750 -0.5976943  0.74769428 0.9899903
IV-III -0.7500 -1.4226943 -0.07730572 0.0244899

Com base no resultado do teste de Tukey, pode-se concluir que a área III é significativamente distinta das áreas I, II e IV, pois para estas três diferenças o valor de p-adj foi menor que 0,05. Outro ponto importante que pode-se concluir com este teste a posteriori é que os ventos na área III são mais fortes que os das áreas I, II e IV, basta observar os valores das diferenças e seus intervalos de confiança para os casos onde as diferenças foram significativas.

Análise de Variância com Dois Fatores

Até o presente momento a aplicação da Análise de Variância se resumiu à um universo amostral bastante enxuto. Assim, o que fazer quando a base de dados é mais extensa? Quando o ``n'' amostral é relativamente maior? A resposta é simples, pois o conceito é o mesmo. O caminho é aplicar os passos igualmente visto na teoria. É necessário verificar os pressupostos da ANOVA, montar as hipóteses a serem testadas, aplicar o teste, compreender seus resultados, e quando necessário, aplicar novos testes para visualizar, da melhor maneira possível, os resultados alcançados.

Deste modo, a base de dados a seguir trata de pesagens de peixes coletados em três áreas distintas. Neste primeiro momento, a questão se resume a verificar se há uma diferença de estrutura corporal (biomassa em kg) entre os peixes das três áreas amostradas. Para isso, a saída é aplicar uma Análise de Variância Simples, considerando todos os pressupostos da análise.

Clique no link ao lado para realizar o download da base de dados - Base de Dados

E agora, se necessitassemos responder se existe uma diferença entre a estrutura corporal (biomassa em kg) destes peixes não somente em função das áreas, mas também em função da profundidade onde os mesmos foram capturados? A alternativa nesta questão ainda seria uma Análise de Variância, porém, agora centrada em dois fatores e não mais em um apenas.

Importando os dados para o R

R source
### Análise de Variância com 2 fatores- importando a base de dados do
### arquivo .csv para o R
dados2 <- read.table("dados2.csv", header = TRUE, sep = ";", dec = ",")

Explorando os dados e suas relações

R source
### Análise Exploratória
par(mfrow = c(1, 2))
plot(dados2$area, dados2$peso, xlab = "Área", ylab = "Peso (Kg)")
plot(dados2$prof, dados2$peso, xlab = "Profundidade", ylab = "Peso (Kg)")
R source
par(mfrow = c(1, 1))

Verificando os pressupostos

  • Pressuposto de normalidade dos dados...

Neste caso, apenas a área "C" não apresentou normalidade dos dados (p-valor < 0,05). Embora o pressuposto de normalidade não tenha sido completamente satisfeito, é conhecida a particularidade da implicação de não normalidade dos dados na ANOVA, sendo assim, decidiu-se manter a aplicação da análise.

R source
### Verificando normalidade dos dados...

### Em função da área...
shapiro.test(dados2$peso[dados2$area == "A"])
R output

    Shapiro-Wilk normality test

data:  dados2$peso[dados2$area == "A"]
W = 0.97687, p-value = 0.3781
R source
shapiro.test(dados2$peso[dados2$area == "B"])
R output

    Shapiro-Wilk normality test

data:  dados2$peso[dados2$area == "B"]
W = 0.96909, p-value = 0.3076
R source
shapiro.test(dados2$peso[dados2$area == "C"])
R output

    Shapiro-Wilk normality test

data:  dados2$peso[dados2$area == "C"]
W = 0.95288, p-value = 0.02472
R source
## Em função da profundidade...
shapiro.test(dados2$peso[dados2$prof == "fundo"])
R output

    Shapiro-Wilk normality test

data:  dados2$peso[dados2$prof == "fundo"]
W = 0.97992, p-value = 0.2634
R source
shapiro.test(dados2$peso[dados2$prof == "raso"])
R output

    Shapiro-Wilk normality test

data:  dados2$peso[dados2$prof == "raso"]
W = 0.98938, p-value = 0.7764

Outra forma de avaliar a normalidade de uma variável é utilizando um Q-Q plot. Com este tipo de avaliação gráfica, é possível identificar onde estão ocorrendo os desvios de normalidade dos dados. Nestes gráficos, pode-se notar que a transgressão da Área C ao pressuposto de normalidade não é tão grave, possibilitando a aplicação do teste.

R source
par(mfrow = c(1, 3))
### Area A
qqnorm(dados2$peso[dados2$area == "A"], main = "Área A")
qqline(dados2$peso[dados2$area == "A"])
### Area A
qqnorm(dados2$peso[dados2$area == "B"], main = "Área B")
qqline(dados2$peso[dados2$area == "B"])
### Area A
qqnorm(dados2$peso[dados2$area == "C"], main = "Área C")
qqline(dados2$peso[dados2$area == "C"])
R source
par(mfrow = c(1, 2))
### Profundidade fundo
qqnorm(dados2$peso[dados2$prof == "fundo"],
       main = "Profundidade - Fundo")
qqline(dados2$peso[dados2$prof == "fundo"])
### Profundidade raso
qqnorm(dados2$peso[dados2$prof == "raso"],
       main = "Profundidade - Raso")
qqline(dados2$peso[dados2$prof == "raso"])
R source
par(mfrow = c(1, 1))
  • Pressuposto de homogeneidade das variâncias...

Com base no teste de Bartlett foi possível perceber que as variâncias, em ambos os tratamentos, são constantes (p-valor > 0,05). Sendo assim, o pressuposto foi atendido, pode-se partir para ANOVA.

R source
### Teste de Homocedasticidade das Variâncias - Teste de Bartlett
bartlett.test(dados2$peso, dados2$area)
R output

    Bartlett test of homogeneity of variances

data:  dados2$peso and dados2$area
Bartlett's K-squared = 1.9953, df = 2, p-value = 0.3687
R source
bartlett.test(dados2$peso, dados2$prof)
R output

    Bartlett test of homogeneity of variances

data:  dados2$peso and dados2$prof
Bartlett's K-squared = 0.31713, df = 1, p-value = 0.5733

Aplicando ANOVA com interação dos fatores

A modelo da ANOVA com interação no R pode ser descrito de duas formas:

  • y ~ x1 * x2 - considerando y a variável dependente em função das xi variáveis explicativas do modelo e a interação dos fatores dada pela multiplicação das variáveis explicativas;

  • y ~ x1 + x2 + x1 : x2 - exatamente o mesmo.

O resultado do modelo 1, com interação dos fatores Área e Profundidade, mostra que a interação dos fatores não é significativa (p-valor > 0,05). Assim, é necessário rodar um novo modelo com os fatores sem interação para que seja possível verificar os efeitos de ambos, independentemente, sobre a variável resposta da análise.

R source
### Aplicando a ANOVA - Modelo linear com interação
mod1 <- lm(dados2$peso ~ dados2$area * dados2$prof)
anova(mod1)
R output
Analysis of Variance Table

Response: dados2$peso
                         Df Sum Sq Mean Sq F value    Pr(>F)    
dados2$area               2 168.12  84.060 15.1203 1.059e-06 ***
dados2$prof               1  13.68  13.681  2.4608    0.1189    
dados2$area:dados2$prof   2  21.76  10.878  1.9567    0.1450    
Residuals               148 822.79   5.559                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neste segundo passo, pode-se observar o efeito significativo do fator área sobre a distribuição dos pesos, no entanto, não há um efeito significativo do fator Profundidade neste caso.

R source
mod2 <- lm(dados2$peso ~ dados2$area + dados2$prof)
anova(mod2)
R output
Analysis of Variance Table

Response: dados2$peso
             Df Sum Sq Mean Sq F value    Pr(>F)    
dados2$area   2 168.12  84.060 14.9299 1.221e-06 ***
dados2$prof   1  13.68  13.681  2.4298    0.1212    
Residuals   150 844.54   5.630                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Por fim, para identificar qual ou quais áreas são distintas entre sí, assim como visto anteriormente, utiliza-se o Teste de Tukey. Com base neste, pode-se perceber que a área A é distinta das áreas B e C, sendo que os pesos nesta é relativamente menor que os pesos nas áreas B e C.

R source
### Teste posterior - Post-Hoc - Teste de Tukey
TukeyHSD(aov(dados2$peso ~ dados2$area))
R output
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = dados2$peso ~ dados2$area)

$`dados2$area`
         diff       lwr      upr     p adj
B-A  2.213228  1.052222 3.374233 0.0000380
C-A  2.172095  1.104962 3.239227 0.0000104
C-B -0.041133 -1.184489 1.102223 0.9960101

Uma alternativa gráfica para visualização das diferenças observadas no teste de Tukey é a aplicação da função plot() ao objeto que armazenou a resultado do teste de Tukey ou diretamente ao teste.

R source
### Visualizando graficamente as diferenças observadas pelo teste de
### Tukey
plot(TukeyHSD(aov(dados2$peso ~ dados2$area)))

Referências

  • Becker, J. L. Estatística Básica: transformando dados em informação. Porto Alegre: Editora Bookman, 488 p. 2015.

  • Casella, G. & Berger, R. L. Inferência Estatística. São Paulo: Cengage Learning, 588 p. 2010.

  • Vieira, S. Análise de Variância (ANOVA). São Paulo: Editora Atlâs, 204 p. 2006.