Capítulo 10 Análise de Variância

Os procedimentos t de duas amostras comparam as médias de duas populações ou as respostas médias a dois tratamentos em um experimento.

Naturalmente, os estudos nem sempre comparam apenas dois grupos. Precisamos de um método que compare qualquer quantidade de médias.

10.1 Comparação de várias médias

Os métodos estatísticos para lidar com comparações múltiplas geralmente apresentam dois passos:

1.Um teste geral para verificarmos se há boa evidência de quaisquer diferenças entre os parâmetros que desejamos comparar.

2.Uma análise de acompanhamento detalhada para decidirmos quais parâmetros são diferentes e estimarmos o tamanho das diferenças.

10.2 A ideia da Análise de Variância

Considere a população com 150 árvores pertencente a um reflorestamento de Mogno Africano (Capítulo 6).

Suponha que retiremos 5 amostras aleatórias desta população com 8 indivíduos em cada amostra.

am1	am2	am3	am4	am5
14.69	15.31	12.91	12.88	10.98
11.96	11.78	9.86	14.57	11.46
11.13	14.55	12.64	14.80	12.49
11.05	14.97	12.79	14.14	15.60
14.24	13.77	13.80	11.80	10.84
16.19	15.64	15.57	12.52	10.96
12.75	7.92	13.88	12.50	13.48
14.15	13.09	16.71	12.09	12.60

Estes dados podem ser encontrados no arquivo dap5.csv

Uma Análise de Variância pode ser executada com as funções lm e anova, como a seguir:

aov_dap5 <- lm(dap ~ amostra, data = dap5)
anova(aov_dap5)

## Analysis of Variance Table
## 
## Response: dap
##           Df  Sum Sq Mean Sq F value Pr(>F)
## amostra    4   7.371  1.8428  0.5075 0.7305
## Residuals 35 127.091  3.6312

Queremos testar a hipótese nula de que não há diferenças entre os diâmetros médios das cinco populações de onde as amostras foram retiradas:

H₀:μ₁ = μ₂ = μ₃ = μ₄ = μ₅

A hipótese alternativa é a de que há alguma diferença, isto é, nem todas as três médias populacionais são iguais:

H₁: nem todas as médias μ são iguais

A priori, sabemos que as cinco amostras foram aleatoriamente retiradas da mesma população. Logo, o teste F da Análise de Variância é não significativo (como esperado). Seu p-valor foi de 0.73 e, por isso, não rejeitamos a hipótese H₀, ou seja, não há evidências que os diâmetros médios das amostras sejam diferentes ou que venham de populações diferentes.

Mas o que ocorreria se houvesse um efeito aditivo em cada uma das amostras? Vamos supor que a amostra 1 tenha um efeito aditivo de +5 unidades no DAP. A amostra 2 de -5, a amostra 3 de +3, a amostra 4 de +2 e, por fim, a amostra 5 não tenha nenhum efeito aditivo. O resultado é este mostrado no arquivo dap5p.csv

am1	am2	am3	am4	am5
19.69	10.31	15.91	14.88	10.98
16.96	6.78	12.86	16.57	11.46
16.13	9.55	15.64	16.80	12.49
16.05	9.97	15.79	16.14	15.60
19.24	8.77	16.80	13.80	10.84
21.19	10.64	18.57	14.52	10.96
17.75	2.92	16.88	14.50	13.48
19.15	8.09	19.71	14.09	12.60

A Análise de Variância, com estes valores atualizados, fica assim:

aov_dap5p <- lm(dap ~ amostra, data = dap5p)
anova(aov_dap5p)

## Analysis of Variance Table
## 
## Response: dap
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## amostra    4 482.71 120.678  33.234 1.791e-11 ***
## Residuals 35 127.09   3.631                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para as mesmas hipóteses anteriores, o teste F da Análise de Variância é significativo (p-valor ~ 0). Neste caso, rejeitamos a hipótese H₀ e podemos concluir que há evidências de que as amostras tenham médias de diâmetro diferentes ou que seja oriundas de populações diferentes.

10.3 Análise de Variância para um fator

Exemplo 10.1 (Análise de Variância para um fator) Pesquisadores estudaram a relação entre variedades da flor tropical Heliconia, na ilha de Dominica, e as diferentes espécies de beija-flores que fertilizam essas flores.

Acredita-se que os comprimentos das flores e as formas dos bicos dos beija-flores evoluíram juntos e se adaptaram uns aos outros.

Se isso for verdade, as variedades de flores fertilizadas por diferentes espécies de beija-flores devem ter diferentes distribuições de comprimentos.

O arquivo bflor.xlsx fornece as medidas de comprimentos (em milímetros, mm) para amostras de três variedades de Heliconia, cada uma fertilizada por uma espécie diferente de beija-flor. As espécies são:

Heliconia bihai
Heliconia caribaea vermelha
Heliconia caribaea amarela

Em particular, os comprimentos médios de suas flores são diferentes?

As três variedades têm distribuições com comprimentos diferentes?

Queremos testar a hipótese nula de que não há diferenças entre os comprimentos médios das três populações de flores:

H₀: μ₁ = μ₂ = μ₃

A hipótese alternativa é a de que há alguma diferença, isto é, nem todas as três médias populacionais são iguais:

H₁: nem todas as μ₁, μ₂ e μ₃ são iguais

Iniciemos com uma análise exploratória:

bflor <- readxl::read_excel("bflor.xlsx") %>%
  mutate(especie = factor(especie))

bflor %>%
  group_by(especie) %>%
  summarise(
    n = n(),
    media = mean(comprimento),
    desvpad = sd(comprimento),
    var = var(comprimento)
  )

## # A tibble: 3 × 5
##   especie     n media desvpad   var
##   <fct>   <int> <dbl>   <dbl> <dbl>
## 1 Hb         16  47.6   1.21  1.47 
## 2 Hca        15  36.2   0.975 0.951
## 3 Hcv        23  39.7   1.80  3.24

boxplot(comprimento ~ especie, data = bflor)

Efetuamos a Análise de Variância propriamente dita:

aov_bflor <- lm(comprimento ~ especie, data = bflor)
anova(aov_bflor)

## Analysis of Variance Table
## 
## Response: comprimento
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## especie    2 1082.87  541.44  259.12 < 2.2e-16 ***
## Residuals 51  106.57    2.09                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CONCLUSÃO: Há forte evidência de que as três variedades de flores não tenham, todas, o mesmo comprimento médio.

O teste F não diz quais das três médias são significantemente diferentes. Aparentemente, pela nossa análise de dados preliminar, as flores da bihai são visivelmente maiores que as da vermelha ou da amarela.

As flores vermelhas e amarelas são muito próximas, mas as vermelhas tendem a ser mais compridas.

10.4 Condições para a ANOVA

Temos k AASs independentes, uma de cada uma das k populações.

Como de costume, o planejamento da produção dos dados é a condição mais importante para inferência. Uma amostragem viesada ou confundimento pode tornar qualquer inferência sem sentido.

Cada uma das k populações tem uma distribuição Normal com média desconhecida

Procedimentos para comparação de médias não são muito sensíveis à falta de Normalidade. A Anova torna-se mais segura à medida que os tamanhos das amostras aumentam. Quando não houver valores atípicos e as distribuições forem aproximadamente simétricas, podemos usar a Anova com segurança.

Todas as populações têm o mesmo desvio-padrão \(\sigma\), de valor desconhecido.

A terceira condição é problemática. Não é fácil verificar a condição de igualdade dos desvios-padrão populacionais. Testes estatísticos de igualdade dos desvios-padrão são tão sensíveis à ausência de Normalidade que, na prática, têm pouco valor.

Mas, qual é a gravidade de os desvios-padrão serem desiguais?

A Anova não é muito sensível a violações da condição, particularmente quando todas as amostras têm tamanhos iguais ou semelhantes e nenhuma delas é muito pequena. Ao planejar um estudo, tente tomar amostras do mesmo tamanho aproximado de todos os grupos que pretende comparar e não utilize amostras muito pequenas.

Certifique-se, antes de fazer a Anova, de que os desvios-padrão amostrais sejam, pelo menos, semelhantes entre si. Como regra prática: maior desvio-padrão não seja o dobro (ou triplo) do menor.

Exemplo 10.2 (Verificação do pressuposto da Normalidade) .

## curtose
aov_bflor %>% residuals() %>% moments::kurtosis()

## [1] 2.349862

## assimetria
aov_bflor %>% residuals() %>% moments::skewness()

## [1] 0.5223018

## Teste de Shapiro-Wilk
aov_bflor %>% residuals() %>% shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.94889, p-value = 0.02227

## Gráfico dos quantis normais
aov_bflor %>% residuals() %>% qqnorm()
aov_bflor %>% residuals() %>% qqline()

## Histograma
aov_bflor %>% residuals() %>% hist()

## Ramo e folhas
aov_bflor %>% residuals() %>% stem()

## 
##   The decimal point is at the |
## 
##   -2 | 3
##   -1 | 98776666553200
##   -0 | 9988877655522111
##    0 | 1355666788999
##    1 | 8
##    2 | 00223557
##    3 | 4

Pela análise do conjunto dos resultados acima, não há evidência de desvio severo da Normalidade.

Exemplo 10.3 (Verificação do pressuposto da homogeneidade das variâncias) .

## razão maior/menor desvio-padrão
bflor %>% group_by(especie) %>% summarise(desvpad=sd(comprimento)) %>%
mutate(razao=max(desvpad)/desvpad)

## # A tibble: 3 × 3
##   especie desvpad razao
##   <fct>     <dbl> <dbl>
## 1 Hb        1.21   1.48
## 2 Hca       0.975  1.84
## 3 Hcv       1.80   1

## boxplot condicional dos resíduos
boxplot(residuals(aov_bflor)~especie, data = bflor)

## resíduos vs ajustados
plot(aov_bflor,1)

## teste de Bartlett
bartlett.test(residuals(aov_bflor)~especie, data=bflor)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(aov_bflor) by especie
## Bartlett's K-squared = 6.4839, df = 2, p-value = 0.03909

## teste de Levene
car::leveneTest(residuals(aov_bflor)~especie, data=bflor)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  2  4.3722 0.01768 *
##       51                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pela análise do conjunto dos resultados acima, não há evidência de não homogeneidade das variâncias.

10.5 Testes de acompanhamento

Para sabermos quais tratamentos diferem entre si, utilizamos um teste de acompanhamento (teste post hoc).

Entre os mais utilizados, serão tratados neste capítulo os seguintes testes:

Teste de Dunnett
Teste de Tukey
Contrastes ortogonais

10.5.1 Teste de Dunnett

Testa os contrastes envolvendo o(s) tratamento(s) testemunha (ou controle ou placebo).

Exemplo 10.4 (Teste de Dunnett) A homeopatia procura utilizar pequenas doses de substâncias altamente diluídas e geralmente perigosas. Cientistas advertem que há pouca evidência que apoie a homeopatia como tratamento eficaz para qualquer condição específica, exceto para uns poucos estudos.

Em um desses estudos, os pesquisadores fizeram e suturaram uma incisão muscular profunda em ratos anestesiados e, então, associaram-nos aleatoriamente a um de cinco tratamentos:

Arnica - dose alta
Arnica - dose baixa
Estafisagria - dose alta
Estafisagria - dose baixa
Placebo

As feridas eram examinadas diariamente para se determinar o tempo (em dias) até sua completa cicatrização. Os valores estão no arquivo homeopatia.xlsx

Há evidência significante de que o tempo de cicatrização dependa do tratamento recebido?

Após efetuada a Análise de Variância e verificados os pressupostos²⁷, segue-se com o teste de Dunnett

dn_remedio <- emmeans::emmeans(aov_hom, ~remedio, contr="dunnett", 
                                ref=which(names(table(hom$remedio))=="Placebo"))


## contrastes
dn_remedio$contrasts

##  contrast                       estimate    SE df t.ratio p.value
##  (Arnica-alta) - Placebo           -5.67 0.442 70 -12.831  <.0001
##  (Arnica-baixa) - Placebo          -5.87 0.442 70 -13.284  <.0001
##  (Estafisagria-alta) - Placebo     -5.73 0.442 70 -12.982  <.0001
##  (Estafisagria-baixa) - Placebo    -5.53 0.442 70 -12.529  <.0001
## 
## P value adjustment: dunnettx method for 4 tests

## intervalo de confiança
dn_remedio$contrasts %>% confint()

##  contrast                       estimate    SE df lower.CL upper.CL
##  (Arnica-alta) - Placebo           -5.67 0.442 70    -6.78    -4.56
##  (Arnica-baixa) - Placebo          -5.87 0.442 70    -6.98    -4.76
##  (Estafisagria-alta) - Placebo     -5.73 0.442 70    -6.84    -4.62
##  (Estafisagria-baixa) - Placebo    -5.53 0.442 70    -6.64    -4.42
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

## plot
dn_remedio$emmeans %>% plot()

Remédio	Tempo de Cicatrização (dias)
Tabela X. Tempo de cicatrização médio (dias) para uma profunda ferida cirúrgica tratada com remédios homeopáticos.
Arnica-baixa	15.47	*
Estafisagria-alta	15.60	*
Arnica-alta	15.67	*
Estafisagria-baixa	15.80	*
Placebo	21.33
Médias seguidas por * diferem do tratamento placebo pelo teste de Dunnett (p<0,05).

Pelo teste de Dunnett, verifica-se que todos os remédios homeopáticos levaram a um tempo de cicatrização menor do que o tratamento placebo.

10.5.2 Teste de Tukey

Criado por John Tukey (1915–2000), realiza todos os contrastes (comparações) possíveis entre as médias populacionais mantendo o nível de significância.

Em lugar do teste t, utiliza o valor crítico m, que depende do número de populações, do número de observações e do nível de significância.

As hipóteses testadas são:

H₀: μ_i = μ_j

H₁: μ_i ≠ μ_j

para todas as médias populacionais

Ainda, há muitos outros teste de comparações múltiplas similares ao teste de Tukey. Se você é capaz de interpretar os resultado deste teste, será capaz de entender os resultados de muitos outros.

Exemplo 10.5 (Teste de Tukey) Voltando ao exemplo anterior das Flores de Heliconia ²⁸.

Utilizaremos o pacote emmeans para o cálculo do teste de Tukey.

tk_especie <- emmeans::emmeans(aov_bflor, ~especie, contr = "tukey")

## contrastes
tk_especie$contrasts

##  contrast  estimate    SE df t.ratio p.value
##  Hb - Hca     11.42 0.520 51  21.977  <.0001
##  Hb - Hcv      7.89 0.471 51  16.759  <.0001
##  Hca - Hcv    -3.53 0.480 51  -7.361  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

## Intervalo de confiança
tk_especie$contrasts %>% confint()

##  contrast  estimate    SE df lower.CL upper.CL
##  Hb - Hca     11.42 0.520 51    10.16    12.67
##  Hb - Hcv      7.89 0.471 51     6.75     9.02
##  Hca - Hcv    -3.53 0.480 51    -4.69    -2.37
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 3 estimates

## letras (CLD)
tk_especie$emmeans %>% multcomp::cld(Letters = letters)

##  especie emmean    SE df lower.CL upper.CL .group
##  Hca       36.2 0.373 51     35.4     36.9  a    
##  Hcv       39.7 0.301 51     39.1     40.3   b   
##  Hb        47.6 0.361 51     46.9     48.3    c  
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

## plot
tk_especie$emmeans %>% plot()

A Tabela abaixo mostra um exemplo de como os resultados do teste podem ser apresentados.

Espécie	Comprimento (mm)
Tabela X. Comprimento médio (mm) de flores de três espécies do gênero Heliconia.
Hca	36.18	a
Hcv	39.71	b
Hb	47.60	c
Médias seguidas por letras iguais não diferem entre si pelo teste de Tukey (p<0,05).

Podemos concluir que as flores da variedade H. bihai tem maior comprimento que as demais. As flores de H. caribaea vermelha são maiores que H. caribaea amarela. Esta última tem o menor comprimento entre todas.

10.5.3 Teste de Scott Knott

É um teste de agrupamento que evita a ambiguidade causada em algumas situações pelo teste de Tukey, principalmente nos estudos com uma grande quantidade de tratamentos.

Utiliza-se o pacote ScottKnott²⁹ para a realização do teste.

Exemplo 10.6 (Teste de Scott Knott) Exemplo de utilização do teste de Scott Knott no exemplo das flores de Heliconia .

sk_bflor <- ScottKnott::SK(aov_bflor) 
sk_bflor$out

## $Result
##     Means G1 G2 G3
## Hb  47.60  a      
## Hcv 39.71     b   
## Hca 36.18        c
## 
## $Sig.level
## [1] 0.05
## 
## $Replicates
## [1] 16 23 15

10.5.4 Contrastes ortogonais

O teste de Tukey (e outros testes de comparações múltiplas) fornecem conclusões sobre todos os contrastes possíveis. Alguns (ou muitos) destes contrastes podem não ser de nosso interesse. Outros contrastes de nosso interesse podem ser deixados de lado.

Se você tiver questões específicas em mente antes de produzir os dados, é muito mais eficiente planejar a análise para responder a estas questões.

Não use os contrastes para testar todas as comparações possíveis!

Contrastes só serão válidos se você tiver alguma razão científica para testar uma hipótese em particular antes mesmo de você coletar os dados.

Exemplo 10.7 (Contrastes ortogonais) Para detectarmos a presença de insetos nocivos nos campos de plantações, podemos colocar cartões cobertos com um material pegajoso e examinar os insetos presos nos cartões.

Quais cores mais atraem os insetos?

Pesquisadores colocaram seis cartões de quatro cores (Amarelo, Azul, Branco e Verde) em locais aleatórios de um campo de aveia e contaram o número de besouros de folhas de cereal apanhados. Os resultados estão no arquivo besouro.xlsx.

Podemos executar o teste de Tukey, mas, antes de coletar os dados, já tínhamos questões específicas em mente:

Suspeitamos que cores quentes (amarelo e verde) são mais atrativas que cores frias (azul e branco)
Branco e azul devem ter respostas semelhantes
Verde e amarelo devem ter respostas semelhantes

Assim, testamos três hipóteses:

Hipótese 1:

H₀: μ_Az = μ_Br

H₁: μ_Az ≠ μ_Br

Hipótese 2:

H₀: μ_Am = μ_Ve

H₁: μ_Am ≠ μ_Ve

Hipótese 3:

H₀:(μ_Am+μ_Ve)=(μ_Az+μ_Br)

H~_1:(μAm_+μVe_)>(μAz_+μBr~)

Duas hipóteses envolvem comparações aos pares. A terceira é mais complexa e unilateral.

Reescrevemos as três hipóteses em termos de três contrastes:

C₁ = (0)(μ_Am) + (+1)(μ_Az) + (−1)(μ_Br) + (0)(μ_Ve)

C₂ = (−1)(μ_Am) + (0)(μ_Az) + (0)(μ_Br) + (+1)(μ_Ve)

C₃ = (+1)(μ_Am) + (−1)(μ_Az) + (−1)(μ_Br) + (+1)(μ_Ve)

Verifique que a soma dos coeficientes em cada linha deve ser igual a 0 (zero) (senão não é um contraste).

Para verificar a ortogonalidade entre os contrastes, a soma do produtos dos coeficientes também deve ser igual a 0 (zero).

As hipóteses para estes contrastes são:

Hipótese 1:

H₀: C₁ = 0

H₁: C₁ ≠ 0

Hipótese 2:

H₀: C₂ = 0

H₁: C₂ ≠ 0

Hipótese 3:

H₀: C₃ = 0

H₁: C₃ > 0

#Contrastes
K.bes <- rbind("C1"=c(  0, +1, -1, 0),
            
               "C2"=c( -1,  0,  0, +1),
           
               "C3"=c( +1, -1, -1, +1))

#verificar se contrastes são ortogonais (soma=0)
sum(K.bes[1,]*K.bes[2,])

## [1] 0

sum(K.bes[1,]*K.bes[3,])

## [1] 0

sum(K.bes[2,]*K.bes[3,])

## [1] 0

library(multcomp)
ctr_besouro <- glht(aov_bes, linfct=mcp(cor=K.bes))
summary(ctr_besouro)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = numero ~ cor, data = bes)
## 
## Linear Hypotheses:
##         Estimate Std. Error t value Pr(>|t|)    
## C1 == 0   -1.333      3.274  -0.407 0.967629    
## C2 == 0  -16.000      3.274  -4.886 0.000255 ***
## C3 == 0   47.333      4.631  10.221  < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

## Intervalo de confiança
ctr_besouro %>% confint()

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = numero ~ cor, data = bes)
## 
## Quantile = 2.5958
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##         Estimate lwr      upr     
## C1 == 0  -1.3333  -9.8332   7.1666
## C2 == 0 -16.0000 -24.4999  -7.5001
## C3 == 0  47.3333  35.3127  59.3540

Estamos 95% confiantes de que o número médio de besouros atraídos pelas cores Verde e Amarelo é superior ao número de besouros atraídos pelas cores Branco e Azul.

Há grande evidência que o contraste C3 é maior que zero.

Os outros contrastes mostram que:

Azul e Branco não diferem entre si
Verde e Amarelo diferem entre si

9 Testes para duas amostras

11 Análise de Variância para 2 fatores