Análise de Variância - Prática e Comparações Múltiplas

EST212 - Bioestatística

Helgem de Souza

Introdução

Antes de iniciar nossa aula:

  1. Abra o RStudio.

  2. Crie um novo script.

  3. Crie uma pasta na área de trabalho com seu nome.

  4. Defina a pasta como diretório de trabalho

  5. Salve o script criado na pasta com o nome "aula12_est212.R".

  6. Baixe do Moodle os arquivos crescimento_plantas.csv e penguins.csv salve na pasta criada.

  7. Leia o arquivo crescimento_plantas.csv para o objeto plantas e o arquivo penguins.csv para o objeto penguins.

Observação: Ao ler os arquivos, na função read.csv, utilize o parâmetro (stringsAsFactors = TRUE). Ele garante que as variáveis de texto sejam importadas como fatores.

Introdução

Na última aula, exploramos os fundamentos teóricos básicos da ANOVA. Nessa aula, nosso foco principal será a prática.

Conforme vimos, a ANOVA tem como objetivo a comparação de uma variável contínua em grupos definidos por uma variável discreta (fator). Ela apresenta as seguintes hipóteses:

  • \(H_0\): As médias de todos os grupos são iguais

  • \(H_1\): Pelo menos um dos grupos apresenta média diferente dos demais.

Introdução

Vimos também alguns requisitos para sua utilização:

  • A variável resposta deve ser numérica, enquanto a variável dependente deve ser qualitativa (controlada no experimento).

  • As observações devem ser independentes entre si e entre os grupos (controlada no experimento).

  • Os resíduos devem apresentar normalidade, salvo o caso de grandes amostras (grupos com mais de 30 elementos)

  • As variâncias dos grupos devem ser estatisticamente iguais, ou seja, deve ocorrer a homogeneidade de variâncias.

Hoje vamos aprender a verificar os requisitos e a executar o teste propriamente dito.

Homogeneidade de variâncias

A homogeneidade de variâncias entre grupos é verificada por meio de um teste semelhante ao teste F, o Teste de Levene.

O Teste de Levene é um teste estatístico cujo objetivo é verificar se dois ou mais grupos apresentam variâncias iguais entre si, ou se existem diferenças significativas. Suas hipóteses são:

  • \(H_0\): Todos os grupos apresentam variâncias iguais

  • \(H_1\): Existe pelo menos um grupo com variância diferente das demais.

Não entraremos em detalhes técnicos sobre a construção do teste. Assim como os demais testes, um p-valor significativo (menor que o nível \(\alpha\) de significância) indica que existem variâncias diferentes entre os grupos.

Teste de Levene

Para executar o teste de Levene no R, utilizamos a seguinte função:

  • leveneTest(formula, data = dados),

em que o objeto formula é a fórmula usual do tipo x ~ ye dados representa o conjunto de dados em estudo.

A função leveneTest pertence ao pacote car. Ou seja, antes de utilizarmos o teste, precisamos carregar o pacote utilizando a função library("car").

Independente do número de testes realizados, precisamos carregar o pacote apenas uma vez.

Para exemplificar sua aplicação, vamos continuar nosso exemplo de crescimento de plantas, iniciado na aula anterior.

Teste de Levene - Exemplo

Na aula anterior, começamos a explorar a relação entre o peso de plantas e o tratamento aplicado, com base no conjunto de dados plantas. Vejamos novamente seu boxplot:

#Boxplot - Peso de plantas por tratamento
boxplot(peso ~ grupo, data = plantas)

Vimos que aparentemente os tratamentos influenciam no peso das plantas (trt1 < ctrl < trt2) e que a variabilidade entre os grupos é visualmente semelhante. Agora vamos aplicar o teste de Levene.

Teste de Levene - Exemplo

Vamos testar ao nível de 5% de significância as seguintes hipóteses:

  • \(H_0\): Todos os grupos apresentam as variâncias de peso iguais.

  • \(H_1\): Existe pelo menos um grupo com variância diferente das demais.

Vejamos como executar o teste

#Carregando o pacote car (necessário fazer apenas uma vez)
library("car")

#Execução do Teste de Levene
leveneTest(peso ~ grupo, data = plantas)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.1192 0.3412
      27

O p-valor obtido foi 0,3412, que é maior que o nível de significância adotado. Logo, temos evidências ao nível de 5% de significância para afirmar que todos os grupos apresentam variâncias iguais.

ANOVA - Aplicação do teste

Agora que sabemos verificar se os grupos apresentam variâncias homogêneas, falta apenas verificar a normalidade dos resíduos. Entretanto, essa verificação é realizada após o execução da ANOVA. Então vamos aprender a realizar o teste e posteriormente a verificar a normalidade dos resíduos.

A aplicação da análise de variância se dá por meio da função

aov(formula, data = dados), em que:

  • formula - novamente nossa fórmula habitual, do tipo x~y

  • dados - representa o conjunto de dados em estudo.

Para acessar os detalhes da análise de variância, armazenaremos seu resultado em um objeto.

Para acessar seus detalhes, utilizaremos novamente a função summary().

Vamos seguir com a análise do nosso exemplo de crescimento de plantas

ANOVA - Aplicação do teste

Agora que sabemos que a variância entre os grupos é homogênea, queremos testar as seguintes hipóteses:

  • \(H_0\): O peso médio das plantas é o mesmo entre os tratamentos.

  • \(H_1\): Existe pelo menos um tratamento cujo peso médio difere dos demais.

Vamos adotar um nível de 5% de significância para o teste. O teste é realizado da seguinte forma:

#Análise de variância - peso por tratamento
anova_peso <- aov(peso ~grupo, data = plantas) #Análise salva no objeto aov_peso

#Resultado da análise de variância
summary(anova_peso)
            Df Sum Sq Mean Sq F value Pr(>F)  
grupo        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Antes de concluirmos, vamos analisar o que significa cada um dos elementos da tabela acima

ANOVA - Aplicação do teste

Analysis of Variance Model
  Df Sum Sq Mean Sq F value Pr(>F)
grupo 2 3.766 1.883 4.846 0.01591
Residuals 27 10.49 0.3886 NA NA

  • Df - graus de liberdade

  • Sum Sq - somas de quadrados (numerador das variâncias entre e dentro)

  • Mean Sq - variância entre e variância dentro, respectivamente

  • F value - valor calculado do teste F

  • Pr(>F) - p-valor do teste

ANOVA - Aplicação do teste

Agora vamos analisar nossas hipóteses, com \(\alpha = 0,05\):

  • \(H_0\): O peso médio das plantas é o mesmo entre os tratamentos.

  • \(H_1\): Existe pelo menos um tratamento cujo peso médio difere dos demais.

# Análise de variância - Peso por grupo
summary(anova_peso)
            Df Sum Sq Mean Sq F value Pr(>F)  
grupo        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como o p-valor é menor que o nível de significância (0,0159 < 0,05), temos evidências, ao nível de 5% de significância que nos levam a não aceitar a hipótese nula. Ou seja, podemos afirmar que pelo menos um dos grupos apresenta média diferente dos demais.

Verificação de normalidade

Agora que executamos o teste, para que todos os requisitos sejam validados e o teste tenha validade estatística, precisamos testar a normalidade dos resíduos.

Os resultados da ANOVA foram armazenados no objeto aov_peso. Os resíduos estão armazenados como variável nesse objeto e podem ser acessados inserindo o sufixo $residuals após o objeto anova_peso:

#Acesso aos resíduos do modelo
anova_peso$residuals
     1      2      3      4      5      6      7      8      9     10     11 
-0.862  0.548  0.148  1.078 -0.532 -0.422  0.138 -0.502  0.298  0.108  0.149 
    12     13     14     15     16     17     18     19     20     21     22 
-0.491 -0.251 -1.071  1.209 -0.831  1.369  0.229 -0.341  0.029  0.784 -0.406 
    23     24     25     26     27     28     29     30 
 0.014 -0.026 -0.156 -0.236 -0.606  0.624  0.274 -0.266

Lembrando que o resíduo é a diferença entre cada observação da variável e a média de seu grupo, ou seja:

\[ e_i = x_{ij} - \mu_j \]

Agora que temos acesso aos resíduos, podemos testar sua normalidade por meio do teste de Shapiro-Wilk, com a função shapiro.test.

Verificação de normalidade

Vamos verificar a normalidade dos resíduos ao nível de 0,1% de significância.

#Teste de Shapiro-Wilk dos resíduos da ANOVA
shapiro.test(anova_peso$residuals)

    Shapiro-Wilk normality test

data:  anova_peso$residuals
W = 0.96607, p-value = 0.4379

Como o p-valor é maior que 0,1, temos evidências que nos levam a não rejeitar a hipótese de normalidade. Portanto, a análise de variância realizada cumpriu todos os requisitos e seus resultados podem ser interpretados com validade estatística.

Em tese, apenas após essa etapa é válido concluir sob rejeitar ou não rejeitar a hipótese nula, apesar de termos feito a interpretação antes dessa verificação.

Caso a normalidade não fosse verificada, as conclusões não seriam válidadas e outro teste, adequado para dados não normais, deveria ser aplicado.

IMPORTANTE: O teste de normalidade pode ser dispensado se todos os grupos apresentarem mais de 30 observações.

Comparações múltiplas

Acabamos de concluir de maneira estatisticamente embasada que o peso médio das plantas é diferente em pelo menos um dos tratamentos.

Essa análise não parece incompleta?

Que pergunta fica em suspenso após essa conclusão?

Comparações múltiplas

Acabamos de concluir de maneira estatisticamente embasada que o peso médio das plantas é diferente em pelo menos um dos tratamentos.

Essa análise não parece incompleta?

Que pergunta fica em suspenso após essa conclusão?

“Se existe pelo menos um grupo com um peso médio diferente, quais seriam essas diferenças?”

Para responder essa pergunta, utilizamos um procedimento de comparações múltiplas.

Um procedimento de comparação múltipla executa testes de comparação 2 a 2 para verificar as eventuais diferenças.

Vamos utilizar o teste de comparações de Tukey

Teste de Comparações de Tukey

O Teste de Tukey é utilizado em comparações 2 a 2. Ou seja, ele irá comparar as médias de todas as possíveis combinações de fatores, 2 a 2. Por exemplo:

  • Se tivermos 3 fatores, teremos \(C^3_2 = \dfrac{3!}{2!(3-2)!} = \dfrac{3!}{2!} = 3\) comparações.

  • Se tivermos 4 fatores, teremos \(C^4_2 = \dfrac{4!}{2!(2-2)!} = \dfrac{4!}{2!2!} = 6\) comparações,

e assim sucessivamente.

Teste de Comparações de Tukey

Para cada comparação, o teste apresenta as seguintes hipóteses:

\(H_0\): A diferença entre as médias dos grupos i e j é igual a zero.

\(H_1\): A diferença entre as médias dos grupos i e j é diferente de zero.

Ou seja, ele verifica se os pares de grupos apresentam a mesma média, ou médias distintas.

O teste de Tukey é realizado por meio da função TukeyHSD, aplicada ao resultado de nossa ANOVA. Vamos verificar por meio de um exemplo como aplicá-lo no R.

Teste de Tukey - Exemplo

Vamos agora finalizar nossa Análise de Variância. Ao analisar a relação entre o peso das plantas e o tratamento aplicado, concluimos que pelo menos uma das médias era diferente das demais.

O resultado da ANOVA foi armazenado no objeto anova_peso. Vamos utilizá-lo na realização do teste de Tukey:

#Teste de Tukey para o peso de plantas
TukeyHSD(anova_peso)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = peso ~ grupo, data = plantas)

$grupo
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

Novamente, antes de interpretar, vamos entender o que cada informação representa:

Teste de Tukey - Exemplo

  • grupo:

      diff lwr upr p adj
    trt1-ctrl -0.371 -1.062 0.3202 0.3909
    trt2-ctrl 0.494 -0.1972 1.185 0.198
    trt2-trt1 0.865 0.1738 1.556 0.01201

  • Cada linha representa uma diferença entre médias. Por exemplo, a linha 1 representa a diferença entre as médias de peso de plantas do grupo trt1 e plantas do grupo crtl.

  • diff representa a diferença entre as médias dos grupos.

  • lwr e upr representam os limites do intervalo de confiança para a diferença de médias.

  • p adj representa o p-valor do teste de hipóteses para a diferença entre os dois grupos.

Teste de Tukey - Exemplo

Agora vamos interpretar os resultados, ao nível de 5% de significância:

#Teste de Tukey para o peso de plantas
TukeyHSD(anova_peso)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = peso ~ grupo, data = plantas)

$grupo
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

  • Não existe diferença significativa entre o tratamento 1 e o grupo controle.

  • Não existe diferença significativa entre o tratamento 2 e o grupo controle

  • Existe diferença significativa entre os tratamentos 1 e 2.

Podemos concluir que plantas tratadas com o tratamento 2 apresentam maior peso, se comparadas com aquelas tratadas com o tratamento 1.

Teste de Tukey - Exemplo

Podemos verificar essas diferenças visualmente, utilizando a função plot:

#Gráfico de comparações múltiplas
plot(TukeyHSD(anova_peso))

Quando o intervalo de confiança das diferenças não contém o zero, temos evidências visuais de que a diferença entre médias é diferente de zero.

Assim, finalizamos todo o processo de análise de variância, desde a identificação da existência ou não de diferenças, até a identificação das diferenças em si.

ANOVA - Passo a passo

Conforme vimos até aqui, para executar uma Análise de Variância, devemos executar os seguintes passos:

  1. Verificar se a variável resposta é contínua e segmentada por uma variável discreta (fator).

  2. Analisar visualmente os dados por meio de estatísticas descritivas como média e desvio padrão e gráficos como o boxplot.

  3. Verificar se a variância entre os grupos é homogênea por meio do Teste de Levene (levene.test).

  4. Executar a ANOVA (aov)

  5. Verificar a normalidade dos resíduos (shapiro.test). Se cada grupo da amostra tiver mais de 30 observações, esse teste é dispensável.

  6. Interpretar o p-valor e concluir se as médias dos grupos são iguais, ou se pelo menos uma delas é diferente das demais.

  7. Caso existam diferenças, executar o teste de Tukey para verificar quais são elas.

Vamos fazer um exercício prático para fixar os conceitos aprendidos até aqui.

Exercício Prático

Vamos utilizar o banco de dados penguins. Esse banco de dados conta com 8 variáveis:

  • espécie: fator com os níveis Adelie, Chinstrap e Gentoo

  • ilha: fator com os níveis Biscoe, Dream e Torgersen

  • comprimento_bico: numérica

  • altura_bico: numérica

  • comprimento_nadadeira: numérica

  • peso: numérica

  • sexo: numérica

  • ano: numérica

Exercício Prático

Queremos verificar se existe diferenças entre o comprimento da nadadeira conforme as espécies:

  • Adelie

  • Chinstrap

  • Gentoo

Queremos testar as seguintes hipóteses:

  • \(H_0\): Os comprimentos médios de nadadeira são iguais entre as espécies

  • \(H_1\): Pelo menos uma das espécies apresenta comprimento médio de nadadeira diferente das demais.

Vamos executar o passo a passo da ANOVA utilizando estes dados.

Passo 1 - Verificação do tipo de variável

Primeiramente precisamos verificar se a variável resposta é numérica e a variável explicativa é discreta (fator). Podemos verificar o tipo de todas as variáveis simultaneamente, bem como seus primeiros valores, utilizando a função str no banco de dados em estudo:

#Tipo de variável
str(penguins)
'data.frame':   344 obs. of  8 variables:
 $ especie              : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ ilha                 : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ comprimento_bico     : num  39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
 $ altura_bico          : num  18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
 $ comprimento_nadadeira: int  181 186 195 NA 193 190 181 195 193 190 ...
 $ peso                 : int  3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
 $ sexo                 : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
 $ ano                  : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

É fácil verificar que a variável comprimento_nadadeira é numérica e a variável espécie é um fator. Logo, os dados atendem ao primeiro requisito.

Passo 2 - Análise Descritiva

Vamos verificar o comportamento do comprimento de nadadeira por espécie, por meio da análise descritiva. Primeiramente verificaremos seu comportamento em termos da média e do desvio padrão, usando a função aggregate:

#Média por grupo
aggregate(comprimento_nadadeira ~ especie, data = penguins, FUN = mean)
    especie comprimento_nadadeira
1    Adelie              189.9536
2 Chinstrap              195.8235
3    Gentoo              217.1870
#Desvio padrão por grupo
aggregate(comprimento_nadadeira ~ especie, data = penguins, FUN = sd)
    especie comprimento_nadadeira
1    Adelie              6.539457
2 Chinstrap              7.131894
3    Gentoo              6.484976

Aparentemente, a espécie Chinstrap possui um comprimento médio maior se comparado às demais, que apresentam comprimento semelhante.

Em relação ao desvio padrão, há indícios de que os grupos apresentam homogeneidade.

Passo 2 - Análise Descritiva

Vamos verificar a relação entre os dados utilizando o boxplot

#Boxplot - Comprimento de nadadeira por espécie
boxplot(comprimento_nadadeira ~ especie, data = penguins,
        main = "Comprimento de nadadeira por espécie", 
        xlab = "Espécie", ylab = "Comprimento de nadadeira")

A relação entre as espécies e o comprimento de nadadeira fica um pouco mais clara ao se observar a análise gráfica, assim como a homogeneidade das variâncias.

Passo 3 - Verificar a homogeneidade de variâncias

A análise descritiva apresenta indícios de que existe homogeneidade de variâncias. Vamos comprovar utilizando o Teste de Levene, ao nível de significância 5% (\(\alpha = 0,05\)). Suas hipóteses são:

  • \(H_0\): Todos os grupos apresentam variâncias iguais

  • \(H_1\): Existe pelo menos um grupo com variância diferente das demais.

#Teste de Levene
leveneTest(comprimento_nadadeira ~ especie, data = penguins)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   2  0.3306 0.7188
      339

Como o p-valor > 0,05, ao nível de 5% de significância, podemos afirmar que as variâncias dos comprimentos de nadadeira entre os grupos são homogêneas.

Passo 4 - Análise de Variância

Como as variâncias são homogêneas, podemos executar a ANOVA. Vamos utilizar 5% de significância. As hipóteses para esse exercício são:

  • \(H_0\): Os comprimentos médios de nadadeira são iguais entre as espécies

  • \(H_1\): Pelo menos uma das espécies apresenta comprimento de nadadeira diferente das demais.

#Análise de Variância
anova_nadadeira <- aov(comprimento_nadadeira ~ especie, data = penguins)

#Detalhes da ANOVA
summary(anova_nadadeira)
             Df Sum Sq Mean Sq F value Pr(>F)    
especie       2  52473   26237   594.8 <2e-16 ***
Residuals   339  14953      44                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 observations deleted due to missingness

Antes de interpretarmos seu resultado, precisamos verificar a normalidade dos resíduos

Passo 5 - Normalidade dos Resíduos

Antes de interpretar os resultados, precisamos verificar a normalidade dos resíduos. Nesse exemplo, como temos mais de 30 indivíduos por grupo, o teste é dispensável. Vamos realizá-lo para fins ilustrativos.

Os resíduos estão armazenados no objeto anova_nadadeira$residuals.

Vamos veriricar a normalidade com o teste de Shapiro-Wilk, ao nível de 10% de significância:

#Teste de normalidade dos resíduos
shapiro.test(anova_nadadeira$residuals)

    Shapiro-Wilk normality test

data:  anova_nadadeira$residuals
W = 0.99452, p-value = 0.2609

Como p-valor > 0,1, temos evidências ao nível de 10% de significância para não rejeitar a hipótese nula, logo, podemos concluir que os resíduos são normais. Portanto, todos os requisitos foram atendidos e agora podemos interpretar o resultado da ANOVA com validade estatística.

Passo 6 - Interpretação do resultado

Vamos novamente examinar o resultado da ANOVA:

#Resultados da ANOVA
summary(anova_nadadeira)
             Df Sum Sq Mean Sq F value Pr(>F)    
especie       2  52473   26237   594.8 <2e-16 ***
Residuals   339  14953      44                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 observations deleted due to missingness

Como p-valor < 0,05, temos evidências ao nível de 5% de significância para não aceitar a hipótese nula. Ou seja, pelo menos uma das espécies apresenta comprimento médio de nadadeira diferente das demais.

Passo 7 - Teste de Tukey

Como identificamos a existência de diferenças nas médias, vamos verificá-las por meio do teste de Tukey, ao nível de 5% de significância

#Teste de Tukey
TukeyHSD(anova_nadadeira)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = comprimento_nadadeira ~ especie, data = penguins)

$especie
                      diff       lwr       upr p adj
Chinstrap-Adelie  5.869887  3.586583  8.153191     0
Gentoo-Adelie    27.233349 25.334376 29.132323     0
Gentoo-Chinstrap 21.363462 19.000841 23.726084     0

Note que todos os p-valores são significativos (menores que 0,05), ou seja, as três espécies de penguins apresentam comprimentos médios de nadadeiras diferentes. Pela análise descritiva, podemos concluir que a ordem crescente de comprimento de nadadeiras é a seguinte: Adelie <Chinstrap < Gentoo.

Exercício

  1. Refaça o exemplo de crescimento de plantas.

  2. Refaça o exercício prático que compara o comprimento da nadadeira por espécie.

  3. Verifique por meio de uma ANOVA se existem diferenças significativas no comprimento de bico dos penguins por espécie.

Não se esqueça de incluir no script a interpretação dos resultados.