Como realizar anova unidirecional em r

Uma ANOVA unidirecional é usada para determinar se há ou não uma diferença estatisticamente significativa entre as médias de três ou mais grupos independentes.

Esse tipo de teste é chamado de ANOVA unidirecional porque analisamos o impacto de uma variável preditora em uma variável de resposta.

Nota : Se, em vez disso, estivéssemos interessados no impacto de duas variáveis preditoras em uma variável de resposta, poderíamos realizar uma ANOVA bidirecional .

Como realizar ANOVA unidirecional em R

O exemplo a seguir ilustra como realizar uma ANOVA unidirecional em R.

Fundo

Suponha que queiramos determinar se três programas de exercícios diferentes têm impactos diferentes na perda de peso. A variável preditora que estudamos é o programa de exercícios e a variável resposta é a perda de peso, medida em libras.

Podemos realizar uma ANOVA unidirecional para determinar se existe uma diferença estatisticamente significativa entre a perda de peso resultante dos três programas.

Recrutamos 90 pessoas para participar de um experimento no qual designamos aleatoriamente 30 pessoas para seguir o Programa A, o Programa B ou o Programa C durante um mês.

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

 #make this example reproducible
set.seed(0)

#create data frame
data <- data.frame(program = rep(c("A", "B", "C"), each = 30),
                   weight_loss = c(runif(30, 0, 3),
                                   runif(30, 0, 5),
                                   runif(30, 1, 7)))

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

# program weight_loss
#1 A 2.6900916
#2 A 0.7965260
#3 A 1.1163717
#4 A 1.7185601
#5 A 2.7246234
#6 A 0.6050458

A primeira coluna do quadro de dados mostra o programa do qual a pessoa participou durante um mês e a segunda coluna mostra a perda total de peso que a pessoa experimentou no final do programa, medida em libras.

Explore os dados

Antes mesmo de ajustar o modelo ANOVA unidirecional, podemos entender melhor os dados encontrando a média e o desvio padrão da perda de peso para cada um dos três programas usando o pacote dplyr :

 #load dplyr package
library (dplyr)

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

# A tibble: 3 x 3
# program mean sd
#      
#1 A 1.58 0.905
#2 B 2.56 1.24 
#3 C 4.13 1.57  

Também podemos criar um boxplot para cada um dos três programas para visualizar a distribuição da perda de peso de cada programa:

 #create boxplots
boxplot(weight_loss ~ program,
data = data,
main = "Weight Loss Distribution by Program",
xlab = "Program",
ylab = "Weight Loss",
col = "steelblue",
border = "black")

A partir destes gráficos de caixa, podemos ver que a perda média de peso é mais elevada para os participantes do Programa C e a perda média de peso é mais baixa para os participantes do Programa A.

Também podemos ver que o desvio padrão (o “comprimento” do boxplot) para perda de peso é um pouco maior no programa C em comparação com os outros dois programas.

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

Ajuste de modelo ANOVA unidirecional

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

aov(variável de resposta ~ preditor_variable, dados = conjunto de dados)

Em nosso exemplo, podemos usar o código a seguir para ajustar o modelo ANOVA unidirecional, usando perda_de_peso como variável de resposta e programa como variável preditora. Podemos então usar a função summary() para exibir o resultado do nosso modelo:

 #fit the one-way ANOVA model
model <- aov(weight_loss ~ program, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#program 2 98.93 49.46 30.83 7.55e-11 ***
#Residuals 87 139.57 1.60                     
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir dos resultados do modelo, podemos perceber que o programa de variáveis preditoras é estatisticamente significativo ao nível de significância de 0,05.

Ou seja, existe uma diferença estatisticamente significativa entre a perda média de peso resultante dos três programas.

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 unidirecional assume:

1. Independência – as observações de cada grupo devem ser independentes umas das outras. Como utilizamos um desenho randomizado (ou seja, designamos aleatoriamente os participantes para os programas de exercícios), essa suposição deve ser atendida para que não tenhamos que nos preocupar muito com isso.

2. Normalidade – a variável dependente deve ter uma distribuição aproximadamente normal para cada nível da variável preditora.

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

Uma maneira de verificar as suposições de normalidade e igualdade de variância é usar a função plot() , que produz quatro gráficos de verificação de modelo. Em particular, estamos particularmente interessados nas duas parcelas a seguir:

• Residuais vs. ajustado – este gráfico mostra a relação entre resíduos e valores ajustados. Podemos usar este gráfico para avaliar aproximadamente se a variação entre os grupos é aproximadamente igual ou não.
• Gráfico QQ – este gráfico exibe os resíduos padronizados em relação aos quantis teóricos. Podemos usar este gráfico para avaliar aproximadamente se a suposição de normalidade é atendida ou não.

O código a seguir pode ser usado para produzir esses gráficos de verificação de modelo:

 plot(model)

O gráfico QQ acima nos permite verificar a suposição de normalidade. Idealmente, os resíduos padronizados estariam ao longo da linha reta diagonal do gráfico. Porém, no gráfico acima podemos observar que os resíduos se desviam um pouco da linha no início e no final. Isso indica que nossa suposição de normalidade pode ser violada.

Os Residuais vs. O gráfico ajustado acima nos permite verificar nossa suposição de variâncias iguais. Idealmente, gostaríamos que os resíduos fossem distribuídos igualmente para cada nível dos valores ajustados.

Podemos ver que os resíduos estão muito mais dispersos para os valores ajustados mais elevados, indicando que a nossa suposição de igualdade de variâncias pode ser violada.

Para testar formalmente variâncias iguais, poderíamos realizar o teste de Levene usando o pacote car :

 #load car package
library (car)

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

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

O valor p do teste é 0,01862 . Se utilizarmos um nível de significância de 0,05, rejeitaríamos a hipótese nula de que as variâncias são iguais entre os três programas. Porém, se utilizarmos um nível de significância de 0,01, não rejeitaremos a hipótese nula.

Embora pudéssemos tentar transformar os dados para garantir que os nossos pressupostos de normalidade e igualdade de variâncias sejam cumpridos, por enquanto não nos preocuparemos muito com isso.

Analise as diferenças de tratamento

Depois de verificarmos que as suposições do modelo são atendidas (ou razoavelmente 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 ~ program, data = data)
#
#$program
# diff lwr upr p adj
#BA 0.9777414 0.1979466 1.757536 0.0100545
#CA 2.5454024 1.7656076 3.325197 0.0000000
#CB 1.5676610 0.7878662 2.347456 0.0000199

O valor p indica se há ou não diferença estatisticamente significativa entre cada programa. Os resultados mostram que existe diferença estatisticamente significativa entre a perda média de peso de cada programa no nível de significância de 0,05.

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

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

Os resultados dos intervalos de confiança são consistentes com os resultados dos testes de hipóteses.

Em particular, podemos observar que nenhum dos intervalos de confiança para a perda média de peso entre os programas contém o valor zero , indicando que existe uma diferença estatisticamente significativa na perda média entre os três programas.

Isso é consistente com todos os valores de p para nossos testes de hipótese sendo inferiores a 0,05.

Relatório de resultados de ANOVA unidirecional

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

Uma ANOVA unidirecional foi realizada para examinar os efeitos do programa de exercícios   na perda de peso (medido em libras). Houve diferença estatisticamente significativa entre os efeitos dos três programas na perda de peso (F(2, 87) = 30,83, p = 7,55e-11). Foram realizados testes post-hoc de Tukey HSD.

A perda média de peso dos participantes do programa C é significativamente maior que a perda média de peso dos participantes do programa B (p < 0,0001).

A perda média de peso dos participantes do programa C é significativamente maior que a perda média de peso dos participantes do programa A (p < 0,0001).

Além disso, a perda média de peso dos participantes do programa B foi significativamente maior do que a perda média de peso dos participantes do programa A (p = 0,01).

Recursos adicionais

Os tutoriais a seguir fornecem informações adicionais sobre ANOVAs unilaterais:

Uma introdução à ANOVA unidirecional
Um guia para usar testes post-hoc com ANOVA
O guia completo: como relatar resultados de ANOVA