Como realizar anova bidirecional em r

Uma ANOVA bidirecional (“análise de variância”) é usada para determinar se há ou não uma diferença estatisticamente significativa entre as médias de três ou mais grupos independentes que foram divididos em dois fatores.

Este tutorial explica como realizar uma ANOVA bidirecional em R.

Exemplo: ANOVA bidirecional em R

Digamos que queremos determinar se a intensidade do exercício e o gênero impactam a perda de peso. Neste caso, os dois factores que estamos a analisar são o exercício e o género , e a variável de resposta é a perda de peso, medida em libras.

Podemos realizar uma ANOVA bidirecional para determinar se o exercício e o gênero impactam a perda de peso e para determinar se existe uma interação entre o exercício e o gênero na perda de peso.

Estamos recrutando 30 homens e 30 mulheres para participar de um experimento no qual designamos aleatoriamente 10 de cada um para seguir um programa sem exercícios, com exercícios leves ou com exercícios intensos por um mês.

O código a seguir cria o quadro de dados com o qual trabalharemos:

 #make this example reproducible
set.seed(10)

#create data frame
data <- data.frame(gender = rep(c("Male", "Female"), each = 30),
                   exercise = rep(c("None", "Light", "Intense"), each = 10, times = 2),
                   weight_loss = c(runif(10, -3, 3), runif(10, 0, 5), runif(10, 5, 9),
                                   runif(10, -4, 2), runif(10, 0, 3), runif(10, 3, 8)))

#view first six rows of data frame
head(data)

# gender exercise weight_loss
#1 Male None 0.04486922
#2 Male None -1.15938896
#3 Male None -0.43855400
#4 Male None 1.15861249
#5 Male None -2.48918419
#6 Male None -1.64738030

#see how many participants are in each group
table(data$gender, data$exercise)

# Intense Light None
# Female 10 10 10
# Male 10 10 10

Explore os dados

Antes mesmo de ajustar o modelo ANOVA bidirecional, podemos entender melhor os dados encontrando a média e o desvio padrão da perda de peso para cada um dos seis grupos de tratamento usando o pacote dplyr :

 #load dplyr package
library(dplyr)

#find mean and standard deviation of weight loss for each treatment group
data %>%
  group_by (gender, exercise) %>%
  summarize (mean = mean(weight_loss),
            sd = sd(weight_loss))

# A tibble: 6 x 4
# Groups: gender [2]
# gender exercise means sd
#          
#1 Female Intense 5.31 1.02 
#2 Female Light 0.920 0.835
#3 Female None -0.501 1.77 
#4 Male Intense 7.37 0.928
#5 Male Light 2.13 1.22 
#6 Male None -0.698 1.12 

Também podemos criar um boxplot para cada um dos seis grupos de tratamento para visualizar a distribuição da perda de peso para cada grupo:

 #set margins so that axis labels on boxplot don't get cut off
by(mar=c(8, 4.1, 4.1, 2.1))

#create boxplots
boxplot(weight_loss ~ gender:exercise,
data = data,
main = "Weight Loss Distribution by Group",
xlab = "Group",
ylab = "Weight Loss",
col = "steelblue",
border = "black", 
las = 2 #make x-axis labels perpendicular
)

Podemos constatar de imediato que os dois grupos que participaram em exercício intenso parecem apresentar valores de perda de peso mais elevados. Podemos verificar também que os homens tendem a ter valores de perda de peso mais elevados do que as mulheres, tanto no grupo de exercício intenso como no de exercício leve .

A seguir, ajustaremos o modelo ANOVA bidirecional aos nossos dados para ver se essas diferenças visuais são realmente estatisticamente significativas.

Ajustando o modelo ANOVA bidirecional

A sintaxe geral para ajustar um modelo ANOVA bidirecional em R é:

aov(variável de resposta ~predictor_variable1 *predictor_variable2, dados = conjunto de dados)

Observe que o * entre as duas variáveis preditoras indica que também queremos testar um efeito de interação entre as duas variáveis preditoras.

Em nosso exemplo, podemos usar o código a seguir para ajustar o modelo ANOVA bidirecional, usando perda_de_peso como variável de resposta e gênero e exercício como as duas variáveis preditoras.

Podemos então usar a função summary() para exibir o resultado do nosso modelo:

 #fit the two-way ANOVA model
model <- aov(weight_loss ~ gender * exercise, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#gender 1 15.8 15.80 11.197 0.0015 ** 
#exercise 2 505.6 252.78 179.087 <2e-16 ***
#gender:exercise 2 13.0 6.51 4.615 0.0141 *  
#Residuals 54 76.2 1.41                   
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir dos resultados do modelo, podemos ver que o género , o exercício e a interação entre as duas variáveis são todos estatisticamente significativos ao nível de significância de 0,05.

Verificando as suposições do modelo

Antes de prosseguirmos, precisamos verificar se as premissas do nosso modelo são atendidas para que os resultados do nosso modelo sejam confiáveis. Em particular, uma ANOVA bidirecional assume:

1. Independência – as observações de cada grupo devem ser independentes umas das outras. Como usamos um desenho aleatório , essa suposição deve ser atendida, portanto não precisamos nos preocupar muito com isso.

2. Normalidade – a variável dependente deve ter distribuição aproximadamente normal para cada combinação de grupos dos dois fatores.

Uma forma de testar esta suposição é criar um histograma dos resíduos do modelo. Se os resíduos forem distribuídos aproximadamente normalmente, esta suposição deve ser satisfeita.

 #define model residuals
reside <- model$residuals

#create histogram of residuals
hist(resid, main = "Histogram of Residuals", xlab = "Residuals", col = "steelblue")

Os resíduos têm distribuição aproximadamente normal, portanto podemos assumir que a suposição de normalidade é atendida.

3. Variância igual – as variâncias de cada grupo são iguais ou aproximadamente iguais.

Uma maneira de verificar essa suposição é realizar um teste de Levene para igualdade de variâncias usando o pacote car :

 #load car package
library(car)

#conduct Levene's Test for equality of variances
leveneTest(weight_loss ~ gender * exercise, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
#group 5 1.8547 0.1177
#54  

Como o valor p do teste é maior que nosso nível de significância de 0,05, podemos assumir que nossa suposição de igualdade de variâncias entre grupos é atendida.

Analise as diferenças de tratamento

Depois de verificarmos que as suposições do modelo foram atendidas, podemos então realizar um teste post hoc para determinar exatamente quais grupos de tratamento diferem uns dos outros.

Para nosso teste post hoc, usaremos a função TukeyHSD() para realizar o teste de Tukey para comparações múltiplas:

 #perform Tukey's Test for multiple comparisons
TukeyHSD(model, conf.level=.95) 

#Tukey multiple comparisons of means
# 95% family-wise confidence level
#
#Fit: aov(formula = weight_loss ~ gender * exercise, data = data)
#
#$gender
# diff lwr upr p adj
#Male-Female 1.026456 0.4114451 1.641467 0.0014967
#
#$exercise
# diff lwr upr p adj
#Light-Intense -4.813064 -5.718493 -3.907635 0.0e+00
#None-Intense -6.938966 -7.844395 -6.033537 0.0e+00
#None-Light -2.125902 -3.031331 -1.220473 1.8e-06
#
#$`gender:exercise`
# diff lwr upr p adj
#Male:Intense-Female:Intense 2.0628297 0.4930588 3.63260067 0.0036746
#Female:Light-Female:Intense -4.3883563 -5.9581272 -2.81858535 0.0000000
#Male:Light-Female:Intense -3.1749419 -4.7447128 -1.60517092 0.0000027
#Female:None-Female:Intense -5.8091131 -7.3788841 -4.23934219 0.0000000
#Male:None-Female:Intense -6.0059891 -7.5757600 -4.43621813 0.0000000
#Female:Light-Male:Intense -6.4511860 -8.0209570 -4.88141508 0.0000000
#Male:Light-Male:Intense -5.2377716 -6.8075425 -3.66800066 0.0000000
#Female:None-Male:Intense -7.8719429 -9.4417138 -6.30217192 0.0000000
#Male:None-Male:Intense -8.0688188 -9.6385897 -6.49904786 0.0000000
#Male:Light-Female:Light 1.2134144 -0.3563565 2.78318536 0.2185439
#Female:None-Female:Light -1.4207568 -2.9905278 0.14901410 0.0974193
#Male:None-Female:Light -1.6176328 -3.1874037 -0.04786184 0.0398106
#Female:None-Male:Light -2.6341713 -4.2039422 -1.06440032 0.0001050
#Male:None-Male:Light -2.8310472 -4.4008181 -1.26127627 0.0000284
#Male:None-Female:None -0.1968759 -1.7666469 1.37289500 0.9990364

O valor p indica se há ou não diferença estatisticamente significativa entre cada grupo.

Por exemplo, na última linha acima, vemos que o grupo de homens sem exercício não apresentou uma diferença estatisticamente significativa na perda de peso em comparação com o grupo de mulheres sem exercício (valor p: 0,990364).

Também podemos visualizar os intervalos de confiança de 95% resultantes do teste de Tukey usando a função plot() em R:

 #set axis margins so labels don't get cut off
by(mar=c(4.1, 13, 4.1, 2.1))

#create confidence interval for each comparison
plot(TukeyHSD(model, conf.level=.95), las = 2)

Relatório de resultados de ANOVA bidirecional

Finalmente, podemos relatar os resultados da ANOVA bidirecional de uma forma que resume os resultados:

Uma ANOVA de dois fatores foi realizada para examinar os efeitos do gênero ( masculino, feminino) e do programa de exercícios (nenhum, leve, intenso) na perda de peso (medido em libras). Houve uma interação estatisticamente significativa entre os efeitos do sexo e do exercício na perda de peso (F(2, 54) = 4,615, p = 0,0141). Foram realizados testes post-hoc de Tukey HSD.

Para os homens, um programa de exercícios intenso resultou em perda de peso significativamente maior do que um programa leve (p < 0,0001) ou nenhum programa de exercícios (p < 0,0001). Além disso, nos homens, uma dieta leve resultou em perda de peso significativamente maior do que nenhum regime de exercícios (p < 0,0001).

Para as mulheres, um programa de exercícios intenso resultou em perda de peso significativamente maior do que um programa leve (p < 0,0001) ou nenhum programa de exercícios (p < 0,0001).

Verificações de normalidade e teste de Levene foram realizadas para verificar se os pressupostos da ANOVA foram atendidos.