Qual melhor análise estatística para fazer comparaçã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]:

As médias dos grupos que não estão sendo comparados são ignoradas
Aumenta a chance de erro tipo I e, portanto, aumenta a probabilidade de observar diferenças significativas por acaso
Quanto mais médias (grupos) existem para comparar, maior é o número de testes estatíticos (e.g. testes t)

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.

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.

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}} \]

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

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

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

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 é 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.

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.

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

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

O two-way ANOVA é uma extensão do one-way ANOVA mas agora adicionamos uma segunda variável independente. Assim como no one-way, o two-way também particiona a variância total em variância sistemática e não-sistemática. Revisando, a variância sistemática é aquela que explica o efeito da variável independente na variável dependente; enquanto a variância não sistemática leva em consideração todas as outras fontes de variação não explicadas pela variável independente. Entretanto, ao contrário do one-way ANOVA, o two-way ANOVA precisa ser analisado de forma diferente considerando os efeitos isolados e combinados das variáveis independentes na variável dependente. As mesmas condições que vimos para o one-way, precisam ser satisfeitas para usar o two-way ANOVA.Este tipo de abordagem tem sido empregada para analisar os efeitos de intervenções terapêuticas por meio dos ensaios clínicos.

Lembre que para o caso de uma única variável independente, usamos um modelo de regressão linear simples com apenas um preditor como vimos antes. Então, para o two-way ANOVA, basta estender aquele modelo para incorporar uma nova variável. De fato, podemos incluir novas variáveis independentes no modelo usando característica aditiva e linea com em um modelo de regressão múltipla.

\[ \hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X_{1} + \hat{\beta{_2}} X_{2} + \epsilon \]

Mas, agora temos de adicionar um termo no modelo para a interação entre as variáveis independentes do modelo.

\[ \hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X_{1} + \hat{\beta{_2}} X_{2} + \hat{\beta{_3}}(X_{1}*X_{2}) + \epsilon \]

Portanto, baseado nos dados e usando um two-way ANOVA, podemos formular três questões a serem respondidas:

Qual é o efeito da variável $X_1$ independente da variável $X_2$?
Qual é o efeito da variável $X_2$ independente da variável $X_1$?
Qual é o efeito combinado ou interação entre as variáveis $X_1$ e $X_2$ ?

Portanto, estas diferentes questões visam investigar efeitos distintos como efeitos principais (de cada variável independente) e efeitos de interação (da combinação das variáveis).

Da mesma forma, três hipóteses nulas vão existir para cada uma das questões acima como:

Não existe diferença nas médias de Y na variável $X_1$
Não existe diferença nas médias de Y na variável de $X_2$
Não existe interação entre $X_1$ e $X_2$ no resultado de Y

No Two-way ANOVA, a variância explicada pelo modelo é consequência não apenas de uma mas sim duas variáveis independentes e da interação entre elas. Portanto, podemos dividir (partilhar) a soma dos quadrados como mostrado abaixo.

Figura. Variâncias separadas no two-way ANOVA: http://www4.uwsp.edu/psych/stat/13/anova-2w.htm

O cálculo do F-ratio é realizado de forma similar ao que foi feito no caso de uma única variável independente (One-way ANOVA). Entretanto agora teremos de calcular três valores de F: um para cada fator e um outro para a interação entre os fatores. Vamos usar o exemplo do livro [1] para aplicar este conceito.

Exemplo. “A researcher wishes to see whether the type of gasoline used and the type of automobile driven have any effect on gasoline consumption. Two types of gasoline, regular and high-octane, will be used, and two types of automobiles, two-wheel- and four-wheeldrive, will be used in each group. There will be two automobiles in each group, for a total of eight automobiles used. Using a two-way analysis of variance, the researcher will perform the following steps.”

Este exemplo ilustra um caso que pode ser usado um ANOVA 2x2, ou seja, 2 fatores (variáveis independentes) com 2 níveis cada resultando em 4 grupos como mostra os dados da Tabela abaixo.

Vamos usar o R para calcular o necessário.

# Groups for 2x2 ANOVA reg2 <- c(26.7,25.2) reg4 <- c(28.6,29.3) high2<- c(32.3,32.8) high4<- c(26.1,24.2) # Grand mean grandM <- mean(c(reg2,reg4,high2,high4)) # Sum of squares # Factor Gas reg <- c(reg2,reg4) high <- c(high2,high4) SSgas <- sum(4*(c(mean(reg)-grandM)^2),4*(c(mean(high)-grandM)^2)) print(paste0("The Sum of squares for the factor type of gasoline is ",toString(round(SSgas,2)))) ## [1] "The Sum of squares for the factor type of gasoline is 3.92" # Factor Automobile type wheel2 <- c(reg2,high2) wheel4 <- c(reg4,high4) SSwheel <- sum(4*(c(mean(wheel2)-grandM)^2),4*(c(mean(wheel4)-grandM)^2)) print(paste0("The Sum of squares for the factor type of automobile is ",toString(round(SSwheel,2)))) ## [1] "The Sum of squares for the factor type of automobile is 9.68" # SSwithin SSw <- sum((reg2-mean(reg2))^2,(reg4-mean(reg4))^2,(high2-mean(high2))^2,(high4-mean(high4))^2) print(paste0("The Sum of squares for the unsystematic variation is ",toString(round(SSw,2)))) ## [1] "The Sum of squares for the unsystematic variation is 3.3" # SStotal SStotal <- sum((reg2-grandM)^2,(reg4-grandM)^2,(high2-grandM)^2,(high4-grandM)^2) print(paste0("The Total Sum of squares is ",toString(round(SStotal,2)))) ## [1] "The Total Sum of squares is 70.98" # SSinteraction SSaxb <- SStotal - SSwheel - SSgas - SSw print(paste0("The Sum of squares for the interaction is ",toString(round(SSaxb,2)))) ## [1] "The Sum of squares for the interaction is 54.08"

Assim como no one-way ANOVA, precisamos calcular os graus de liberdade para cada componente para obter o F-ratio. Como mostra a Tabela abaixo [1].

Onde,

a = número de níveis do fator A
b = número de níveis do fator B
n = número de sujeitos por grupo

# Degrees of Freedom of the numerator a <- 2 b <- 2 n <- 2

O F-ratio é calculado da mesma forma que no one-way ANOVA mas agora temos três valores de F associados a cada um dos fatores e com a interação entre eles. Primeiro calculamos o quadrado médio [3] (mean squares) para cada componente como mostrado abaixo.

\[MS_{A}=\frac{SS_{A}}{a-1}\] \[MS_{B}=\frac{SS_{B}}{b-1}\] \[MS_{AxB}=\frac{SS_{AxB}}{(a-1)(b-1)}\]

# Mean Squares MSgas <- SSgas/(a-1) MSwheel <- SSwheel/(b-1) MSgxw <- SSaxb/((a-1)*(b-1)) MSw <- SSw/(a*b*(n-1)) print(paste0("The MSgas, MSwheel, MSgxw and MSw are respectively ",toString(round(c(MSgas,MSwheel,MSgxw,MSw),3)))) ## [1] "The MSgas, MSwheel, MSgxw and MSw are respectively 3.92, 9.68, 54.08, 0.825"

Perceba que os valores de SS e MS são iguais. Isto se deve porque estamos usando um ANOVA 2x2.

E então, podemos obter o valor de F.

# Degrees of freedom of the denominator dfD <- a*b*(n-1) # as the groups are balanced Fgas <- MSgas/MSw Fgas ## [1] 4.751515 Fwheel <- MSwheel/MSw Fwheel ## [1] 11.73333 Fgxw <- MSgxw/MSw Fgxw ## [1] 65.55152 print(paste0("The Fgas, Fwheel, Fgxw are respectively ",toString(round(c(Fgas,Fwheel,Fgxw),3)))) ## [1] "The Fgas, Fwheel, Fgxw are respectively 4.752, 11.733, 65.552"

Encontrar o valor crítico na distribuição F de acordo com o valor de alpha e os graus de liberdade. Como os graus de liberdade do numerador e do denominador são iguais para os fatores e para a interação, porque é um ANOVA 2x2, podemos calcular apenas um valor crítico.

# Find the critical value alpha = 0.05 cv <- qf(alpha,df1 = 1,df2 = dfD, lower.tail=FALSE) print(paste0("The critical value for F is ",toString(round(cv,3)))) ## [1] "The critical value for F is 7.709"

Portanto, como os valores de Fwheel e Fgxw são maiores que o valor crítico, rejeita-se $H_{0}$ e conclui-se que existe efeito do tipo de automóvel e efeito da interação entre as duas variáveis independentes (gasolina e automóvel) no consumo.

Assim como no one-way ANOVA, podemos usar a função aov no R para resolver diretamente o problema.

# Loading the necessary libraries. Install packages if you haven't done yet, instead. library(psych) library(car) ## ## Attaching package: 'car' ## The following object is masked from 'package:psych': ## ## logit library(lsr) # Making a data frame to input in aov DV <- c(reg2,reg4,high2,high4) gasoline <- c(rep("regular",4),rep("high",4)) automobile <- c(rep("two",2),rep("four",2),rep("two",2),rep("four",2)) data <- data.frame(gasoline,automobile,DV) data ## gasoline automobile DV ## 1 regular two 26.7 ## 2 regular two 25.2 ## 3 regular four 28.6 ## 4 regular four 29.3 ## 5 high two 32.3 ## 6 high two 32.8 ## 7 high four 26.1 ## 8 high four 24.2 # Main effects and interaction analysis aov.data = aov(data$DV ~ data$gasoline*data$automobile) summary(aov.data) ## Df Sum Sq Mean Sq F value Pr(>F) ## data$gasoline 1 3.92 3.92 4.752 0.09477 . ## data$automobile 1 9.68 9.68 11.733 0.02665 * ## data$gasoline:data$automobile 1 54.08 54.08 65.552 0.00126 ** ## Residuals 4 3.30 0.82 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como esperado, os resultados obtidos foram os mesmos de [1].

Para facilitar a interpretação dos resultados, vamos visualizar um plot de interação

# Interaction plot interaction.plot(data$gasoline,data$automobile,data$DV, type = "l",xlab="Type of Gasoline",ylab="Consumption")

Perceba que os níveis do fator tipo de gasolina são mostrados no eixo horizontal e os níveis do fator tipo de automóvel estão representados por diferentes curvas no gráfico. Podemos inverter esta situação como mostrado abaixo.

# Interaction plot interaction.plot(data$automobile,data$gasoline,data$DV, type = "l",xlab="Type of Automobile",ylab="Consumption")

Este tipo de interação é chamada de não-ordinal pois a variável dependente aumenta em um nível do outro fator mas diminui no outro nível. Tipicamente, no caso de uma interação não-ordinal perfeita, a média de ambos os fatores para cada um dos níveis será a mesma e, portanto, não haverá main effects.

Na maioria dos casos, um ANOVA fatorial é conduzido na expectativa de identificar padrões de interação entre as variáveis independentes. Por exemplo, a hipótese de que certas combinações de tratamento serão mais eficazes para alterar a variável dependente. Senão bastava conduzir estudos isolados com apenas uma variável independente. A interpretação do efeito de interação é facilitada conduzindo a análise dos efeitos simples, pois precisamos ir um pouco mais adiante e entender quem está influenciando esta interação. Cada linha no plot de interação representa um efeito simples. Por meio da análise dos efeitos simples é possível detectar quais os níveis das variáveis contribuíram mais para as diferenças observadas.

Interação pode ser definida quando os efeitos de uma variável independente mudam nos diferentes níveis da outra variável independente. Então vamos fazer isto no R.

# Conduct simple effects analysis (of Gasoline at each level of Automobile) AB1 <- subset(data, data$automobile == "two") AB1 ## gasoline automobile DV ## 1 regular two 26.7 ## 2 regular two 25.2 ## 5 high two 32.3 ## 6 high two 32.8 aov.AB1 <- aov(AB1$DV ~ AB1$gasoline) summary(aov.AB1) ## Df Sum Sq Mean Sq F value Pr(>F) ## AB1$gasoline 1 43.56 43.56 69.7 0.014 * ## Residuals 2 1.25 0.62 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 AB2 <- subset(data, data$automobile == "four") AB2 ## gasoline automobile DV ## 3 regular four 28.6 ## 4 regular four 29.3 ## 7 high four 26.1 ## 8 high four 24.2 aov.AB2 <- aov(AB2$DV ~ AB2$gasoline) summary(aov.AB2) ## Df Sum Sq Mean Sq F value Pr(>F) ## AB2$gasoline 1 14.44 14.440 14.09 0.0642 . ## Residuals 2 2.05 1.025 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neste caso, podemos observar que houve efeito do tipo de gasolina no two wheel car mas não no four wheel explicando assim a interação observada. Entretanto, nem sempre a interpretação é simples assim. Em algumas situações, podemos obter efeito simples em ambas as situações e ainda assim haver interação. Isto ocorre pois o valor de F não é adequado para indicar o tamanho do efeito.