Introdução

Até aqui vimos como é possível comparar médias de até dois grupos usando o teste t, por exemplo. Entretanto, em muitas situações estamos interessados em examinar se 3 ou mais grupos (ou condições) difererem entre si. Por exemplo, suponha que três turmas de Bioestatística estão sendo ofertadas neste quadrimestre e estamos interessados em saber se o desempenho dos alunos destas turmas foi diferente? Alguém poderia sugerir o seguinte “por que não usar vários testes t para investigar isto?”. Em outras palavras, se comparararmos as turmas A e B; depois as turmas A e C e, finalmente, as turmas B e C não levaria ao que estamos interessados? Embora estes testes consigam cumprir estes objetivos e testar as diferenças, esta abordagem não é adequada por algumas razões [1]:

Portanto, uma outra abordagem é necessária para responder se estas diferenças existem. Este é um problema frequente em estudos clínicos uma vez que necessitamos entender os efeitos de intervenções (e.g. intervenções terapêuticas ou preventivas) quando mais de 2 grupos ou mais de uma fator são envolvidos. Portanto, a análise de variância (ANOVA) é uma das ferramentas mais utilizadas neste tipo de estudo pois é uma abordagem estatística capaz de examinar diferenças observadas nas médias dos grupos (baseado em suas variâncias esperadas) e as variâncias não explicadas devido a chance por exemplo. Então, ao invés de usarmos a distribuição z ou t usaremos uma outra família de distribuição conhecida como F também conhecida como distribuição de Fisher-Snedecor devido a contribuição de Ronald Fisher.

ANOVA como regressão linear

Na verdade esta abordagem não é diferente da que vimos anteriormente em regressão linear onde podemos predizer a variável dependente a partir de uma variável independente (ou a partir de vários preditores se for uma regressão múltipla). A diferença é que em regressão linear simples, as variáveis independente e dependente são variáveis contínuas e o modelo de regressão é expresso da seguinte forma

\[ \hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X + \epsilon \]

Por exemplo, suponha que queremos predizer o salário dos professores (variável dependente, Y) pelo tempo de trabalho em anos (variável independente, X). Note que ambas são variáveis contínuas. Imagine agora que, ao invés do tempo de trabalho, usássemos o gênero (homem vs. mulher) como variável independente, ou seja, uma variável categórica. A Figura abaixo ilustra esta condição. Suponha que 65 (mulher) e 90 (homem) são os níveis da variável independente categórica. Podemos observar que em cada nível a variável dependente (salário, Y) segue uma distribuição normal com média 0 e variância \(\sigma{^2}\). Em outras palavras, para cada nível da variável independente (X) a variável dependente tem média (\(\hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X\)) e desvio padrão (\(\sigma\)).

Figura. Modelo de regressão destacando a distribuição da variável dependente em cada nível da variável independente. Fonte:http://reliawiki.org/index.php/Simple_Linear_Regression_Analysis

Então, ANOVA é um caso especial de análise de regressão linear. De qualquer forma para o caso acima, o teste t seria adequado pois existem apenas dois níveis da variável independente. O ANOVA é útil quando temos mais do que dois níveis na variável independente.

Teste F

O teste F é usado para comparar variâncias que é a base para a ANOVA que, por sua vez, é a técnica usada para comparar se as médias de 3 ou mais grupos são diferentes. A ideia do teste F é usar o resultado do F-ratio para determinar se rejeita ou não a hipótese nula com base na distribuição F. F-ratio, por sua vez, é a razão da variância explicada, ou seja, aquela que pode ser atribuída pela intervenção pela variância não explicada, ou aquela devido a outros fatores que não diretamente relacionadas à intervenção (e.g. erros experimentais).

\[ F = \frac{\text{Variância explicada}}{\text{Variância não explicada}} \] \[ F = \frac{\text{Variância entre grupos}}{\text{Variância intra-grupo}} \]

Portanto, se a variação devido a uma intervenção foi substancialmente maior que a variação não explicada a tendência é que o valor do F-ratio seja maior que 1. O teste F é sempre unicaudal para a direita (positivo) pois a razão não pode ser negativa. O F-ratio é usado para obter o valor crítico da distribuição F e assim aceitar ou refutar a hipótese nula do estudo.

O numerador na equação representa a variância entre grupos e é calculada a partir da razão entre a soma dos quadrados entre os grupos e respectivo graus de liberdade.

\[ \textit{s}^{2}_{B} = \frac{\sum n_{i}(\bar{X}_{i} - \bar{X}_{GM})^2}{\textit{k}-1} \]

A equação acima tenta capturar a variação entre os grupos pois, como os grupos foram submetidos a condições distintas, é esperado uma variação. Por isto é referido como variância sistemática ou explicada.

Por sua vez, a variância intragrupo é calculada pela razão entre a soma dos quadrados intragrupo e do respectivo graus de liberdade.

\[ \textit{s}^{2}_{W} = \frac{\sum (n_{i}-1) \textit{s}^{2}_{i}}{\sum (n_{i}-1)} \]

E então o F-ratio pode ser calculado,

\[ \textit{F} = \frac{\textit{s}^{2}_{B}}{\textit{s}^{2}_{W}} \]

Graus de liberdade

Os graus de liberdade do teste F tem relação com o numerador e o denominador desta razão (variância entre os grupos e intragrupos). Então os graus de liberdade do numerador é k-1 onde k é o número de grupos. E os graus de liberdade do denominador é N-k, onde N é a amostra total. Assim, para encontrar o valor crítico na distribuição F, é preciso saber os graus de liberdade (numerador e denominador) e o nível de significância \(\alpha\).

Um forma de interpretar a variância intragrupo (que o Andy Conway chama de variância não sistemática) é o desvio de cada sujeito da média do grupo que ele pertence. Por que indivíduos estão variando se eles foram submetidos às mesmas condições (e.g. intervenção)? Não podemos explicar esta variação, por isso consideramos não explicada (não sistemática) ou variações atribuídas a chance.

Estas duas variâncias são também conhecidas como médias quadráticas (mean squares) entre grupos (\(MS_{B}\)) e intragrupo (\(MS_{W}\)).

A excelente ilustração abaixo [2] resume estes conceitos que tratamos até aqui portanto se entendermos a Figura é um passo significativo para entender ANOVA.

Figura. (a) Between- and within-group variance is calculated from SSB, the between treatment sum of squares, and SSW, the within treatment sum of squares.. Deviations are shown as horizontal lines extending from grand and sample means. The test statistic, F, is the ratio mean squares MSB and MSW, which are SSB and SSW weighted by d.f. (b) Distribution of F, which becomes approximately normal as k and N increase, shown for k = 3, 5 and 10 samples each of size n = 6. N = kn is the total number of sample values. (c) ANOVA analysis of sample sets with decreasing within-group variance (σw2 = 6,2,1). MSB = 6 in each case. Error bars, s.d. Fonte: http://www.nature.com/nmeth/index.html

Teste de Hipóteses

Assim como nos testes estatísticos anteriores, a ANOVA é um teste de hipótese e, portanto, é necessário declarar as hipóteses do estudo a priori. Entretanto, diferente do que vimos no teste t, temos mais que dois grupos e, portanto, as hipóteses nula e alternativa são declaradas de forma um pouco diferente.

\[ H_0: \mu_{1} = \mu_{2} = ... \mu_{k}\] \[ H_1: \text{Ao menos uma diferença existe entre os grupos comparados}\]

Condições para usar ANOVA

  • A variável dependente é contínua
  • A variável dependente segue uma distribuição normal
  • Homogeneidade das variâncias

Exemplo One-way ANOVA

A melhor maneira de entender como o ANOVA é calculado é por meio de exemplos. Vamos seguir o exemplo do livro [2].

“A researcher wishes to try three different techniques to lower the blood pressure of individuals diagnosed with high blood pressure. The subjects are randomly assigned to three groups; the first group takes medication, the second group exercises, and the third group follows a special diet. After four weeks, the reduction in each person’s blood pressure is recorded. At \(\alpha=\) 0.05, test the claim that there is no difference among the means. The data are shown.”

Vamos usar o R para resolver.

# Input data
medication<-c(10,12,9,15,13)
exercise<-c(6,8,3,0,2)
diet<-c(5,9,12,8,4)

n1 <-5 # sample size 1
n2 <-5 # sample size 2
n3 <-5 # sample size 3
k <- 3 #number of groups
df1 <- k-1 # dof between group
df2 <- n1+n2+n3-k # dof within group
alpha <- 0.05 # sig. level

# Within group mean
medicationM <- mean(medication)
exerciseM <- mean(exercise)
dietM <- mean(diet)

# Calculating grand mean
grandM <- mean(c(medicationM,exerciseM,dietM))

# Between-group variance
sSqBG <- sum(n1*(medicationM-grandM)^2,n2*(exerciseM-grandM)^2,n3*(dietM-grandM)^2)/(k-1)

# Within-group variance
sSqWG <- sum((n1-1)*var(medication),(n2-1)*var(exercise),(n3-1)*var(diet))/sum(n1-1,n2-1,n3-1)

# Find the F-ratio value
Fratio <- sSqBG/sSqWG

# Find the critical value
cv <- qf(alpha,df1,df2, lower.tail=FALSE)
print(paste0("The critical value is ",toString(round(cv,2))))
## [1] "The critical value is 3.89"
print(paste0("The F test value is ",toString(round(Fratio,2))))
## [1] "The F test value is 9.17"

Como podemos observar o valor do teste F é maior que o valor crítico, portanto podemos rejeitar a hipótese nula. Entretanto, perceba que não podemos dizer ainda onde as diferenças se encontram uma vez existem três grupos. Para conseguir determinar onde as diferenças se encontram é necessário usar uma outra abordagem como veremos adiante.

Ao invés de fazer na mão, podemos usar as bibliotecas e funções disponíveis para isto no R.

# Creating independent group labels
ivLabel <- c(rep("med",5),rep("exe",5),rep("die",5))
dvdata  <- c(medication,exercise,diet)

# Performing One-way ANOVA
lmout <- lm(dvdata~ivLabel)
Fratio2 <- summary(lmout)$fstatistic[1]

print(paste0("The F test value is ",toString(round(Fratio2,2))))
## [1] "The F test value is 9.17"
# Alternatively, one can use
aovOut <- aov(dvdata~ivLabel)
out <- summary(aovOut)
# And obviously the result should be the same
print(paste0("The F test value is ",toString(round(out[[1]]$F[1],2))))
## [1] "The F test value is 9.17"

Portanto o mesmo resultado foi obtido nas três abordagens.

Em geral, os resultados do ANOVA são resumidos em uma tabela como mostrado abaixo [2]:

Vamos ver a Tabela de saída do ANOVA como saída do aov.

summary(aovOut)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## ivLabel      2  160.1   80.07   9.168 0.00383 **
## Residuals   12  104.8    8.73                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tamanho do Efeito

O tamanho do efeito fornece uma medida de quanto a variável dependente foi afetada pela variável independente no caso de um estudo experimental. Por exemplo, em um ensaio clínico gostaria de saber qual o efeito que pode ser atribuído a um tratamento. Em outras palavras, o tamanho do efeito examina quanto da variância na variável dependente foi resultado da variável independente. O tamanho do efeito no ANOVA é medido calculando-se o \(\eta^2\) (eta squared) da seguinte forma. \[\eta^2=\frac{\text{Variância explicada}}{\text{Variância total}}=\frac{SS_{B}}{SS_{W}+SS_{B}}\] Via de regra, na interpretação do tamanho do efeito a seguinte recomendação pode ser adotada: pequeno (\(\eta^2\)=0.01), médio(\(\eta^2\)=0.06) e grande (\(\eta^2\)=0.14). Podemos calcular o \(\eta^2\) para o exemplo anterior usando o output do ANOVA.

# Effect size for ANOVA
ss = summary(aovOut)[[1]]$"Sum Sq"
eta.sq = ss[1]/(ss[1] + ss[2])
print(paste0("The eta-squared is ",toString(round(eta.sq,3))))
## [1] "The eta-squared is 0.604"

Portanto, o tamanho do efeito pode ser considerado grande.

Análise a posteriori (Post hoc analysis)

Análise a posteriori é quando um procedimento estatístico é realizado para analisar padrões nos dados, que não foram determinados a priori, após o experimento ter sido concluído. Portanto, é recomendável não fazer este tipo de procedimento antes de realizar o ANOVA, ou seja, não se faz comparações múltiplas antes, mas depois, caso você tenha obtido um valor significativo no ANOVA. É como se conduzíssemos vários testes t para determinar qual tipo de intervenção (medication, exercise, diet) é mais eficaz para reduzir a pressão arterial ao invés de usar ANOVA. Já havíamos comentado que tal procedimento não é adequado pois, quando se faz múltiplas comparações simples, ignoramos as informações dos grupos deixados de lado nesta comparação. Além disso, a chance de cometermos erro do tipo I cresce a medida que mais comparações (testes) são feitos. Esta condição foi observada por John Tukey e nomeada como experimentwise error rate que é a probabilidade de fazer descobertas falsas quando múltiplas comparações são realizadas. Por outro lado, a análise post hoc é necessária para entender onde as diferenças se encontram caso obtenhamos efeito no ANOVA quando existem mais de dois grupos. Existem diversas abordagens para realizar esta análise a posteriori desde as mais conservadoras que de certa forma levam em consideração a probabilidade de cometer este erro, até as menos conservadoras que não levam em consideração. Vamos ver alguns destes métodos a seguir.

Teste de Scheffé

O teste de Scheffé é similar em forma ao teste t mas ele possui uma correção implementada nele para evitar a inflação do erro tipo I independente de quantos testes forem realizados.

\[ \textit{F}_{S} = \frac{(\bar{X}_i-\bar{X}_j)^2}{(\textit{k}-1)MS_{W}(\frac{1}{n_{i}}+\frac{1}{n_{j}})} \]

# Performing Scheffe test
Fs1 <- (mean(dvdata[1:5])-mean(dvdata[6:10]))^2/((k-1)*sSqWG*((1/5)+(1/5)))
Fs2 <- (mean(dvdata[1:5])-mean(dvdata[11:15]))^2/((k-1)*sSqWG*((1/5)+(1/5)))
Fs3 <- (mean(dvdata[11:15])-mean(dvdata[6:10]))^2/((k-1)*sSqWG*((1/5)+(1/5)))

print(paste0("The Schaffe test (Fs) for MedvsExe, MedvsDie and DievsExe were ",toString(round(c(Fs1,Fs2,Fs3),3))))
## [1] "The Schaffe test (Fs) for MedvsExe, MedvsDie and DievsExe were 9.16, 2.525, 2.067"
# Find the critical value
cv <- qf(alpha,2,12, lower.tail=FALSE)
print(paste0("The critical value is ",toString(round(cv,2))))
## [1] "The critical value is 3.89"

Portanto, apenas a comparação entre o os grupos medication vs. exercise foi estatisticamente significativa.

Teste de Tukey

Assim como o Scheffé, o teste de Tukey também é similar ao teste t mas tenta levar em consideração a inflação do erro tipo I ao fazer comparações múltiplas. Este teste também é conhecido como Tukey Honest Significant Difference (HSD) A fórmula usada para este teste é a seguinte.

\[ \textit{q} = \frac{(\bar{X}_i-\bar{X}_j)}{\sqrt{\frac{MS_{W}^2}{n}}} \]

Podemos usar a função lm do R com o TukeyHSD

# Tukey Honest Significant Test (HSD)
q <- TukeyHSD(aov(dvdata~ivLabel))
q
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = dvdata ~ ivLabel)
## 
## $ivLabel
##         diff        lwr      upr     p adj
## exe-die -3.8 -8.7863602  1.18636 0.1465881
## med-die  4.2 -0.7863602  9.18636 0.1031329
## med-exe  8.0  3.0136398 12.98636 0.0028351

Correção de Bonferroni

O teste de Bonferroni é bastante popular e tem como objetivo ajustar o valor de significância do teste baseado no número de comparações (m) para diminuir a chance de erro tipo I. No exemplo que temos trabalhado \(m=3\) e \(\alpha=0.05\), então rejeitaríamos \(H_0\) somente se \(p \leq \frac{\alpha}{\textit{m}}\). Vamos fazer a mesma comparação mas usando a função do R.

# Bonferroni correction
pairwise.t.test(dvdata,ivLabel,p.adjust.method="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dvdata and ivLabel 
## 
##     die    exe   
## exe 0.1943 -     
## med 0.1327 0.0032
## 
## P value adjustment method: bonferroni

Podemos ainda fazer comparações múltiplas sem levar em consideração a inflação do erro tipo I devido aos múltiplos experimentos. Quando fazemos isso, estamos usando o chamado Fisher Least Significance Difference (LSD). No R podemos usar a mesma função anterior mas sem nenhum ajuste (“none”).

# Fisher Least Significance Difference (LSD)
pairwise.t.test(dvdata,ivLabel,p.adjust.method="none")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dvdata and ivLabel 
## 
##     die    exe   
## exe 0.0648 -     
## med 0.0442 0.0011
## 
## P value adjustment method: none

Perceba que os resultados contrastaram com os anteriores (Scheffe, Tukey e Bonferroni) mostrando que o Fisher’s LSD é o mais liberal de todos e permitiu mostrar diferenças quando medication vs. exercise e medication vs. diet foram comparados. Nos outros testes, foi observado diferenças apenas na comparação medication vs. exercise. Note também que os valores de p do Fisher’s LSD são extamente 3 vezes menores que os obtidos no Bonferroni.

Em conclusão, não existe um teste melhor que o outro sendo possível notar que todos eles tem sido usados na literatura. Em geral, é dito que o Fisher’s LSD é o mais liberal, o Bonferroni é o mais conservador e o Tukey é o mais honesto. Existem outros testes com finalidades similares que não veremos desta vez, veja aqui.

A figura abaixo mostra um exemplo de como aplicar One-way ANOVA para comparar a média de três grupos [2].

Figura. (a) Three samples drawn from normal distributions with swit2 = 2 and treatment means mA = 9, mB = 10 and mC = 11. (b) Depiction of deviations with corresponding SS and MS values. (c) Sample means and their differences. P values for paired sample comparison are adjusted for multiple comparison using Tukey’s method. Error bars, 95% CI. Fonte: http://www.nature.com/nmeth/index.html

Two-way ANOVA

Em algumas situações, estamos interessados em estudar o efeito de mais de uma variável independente como preditora da variável dependente. Neste caso, usamos um outro fator. Veremos na próxima aula.

References

  1. Bluman, Allan G. Elementary statistics : a step by step approach / Allan Bluman. — 8th ed.
  2. Altman & Krzywinski (2015).Points of Significance: Simple linear regression. Nature Methods 12, 999–1000.