如何在 r 中执行双向方差分析

双向方差分析（“方差分析”）用于确定跨两个因素划分的三个或更多独立组的平均值之间是否存在统计显着差异。

本教程介绍如何在 R 中执行双向方差分析。

示例：R 中的双向方差分析

假设我们想要确定运动强度和性别是否会影响减肥。在本例中，我们关注的两个因素是运动和性别，响应变量是体重减轻（以磅为单位）。

我们可以进行双向方差分析来确定运动和性别是否影响减肥，并确定运动和性别对减肥是否存在交互作用。

我们正在招募 30 名男性和 30 名女性参加一项实验，我们随机分配每人 10 人，分别进行一个月的无运动、轻度运动或剧烈运动计划。

以下代码创建我们将使用的数据框：

 #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

探索数据

在拟合双向方差分析模型之前，我们可以通过使用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 

我们还可以为六个治疗组中的每一个创建一个箱线图，以可视化每组的体重减轻分布：

 #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
)

我们立即可以看到，参加剧烈运动的两组似乎都有更高的减肥值。我们还可以看到，无论是剧烈运动组还是轻度运动组，男性的减肥值往往高于女性。

接下来，我们将双向方差分析模型与我们的数据进行拟合，看看这些视觉差异实际上是否具有统计显着性。

拟合双向方差分析模型

在 R 中拟合双向方差分析模型的一般语法是：

aov(响应变量 ~predictor_variable1 *predictor_variable2, data = 数据集)

请注意，两个预测变量之间的*表示我们还想测试两个预测变量之间的交互作用。

在我们的示例中，我们可以使用以下代码来拟合双向方差分析模型，使用weight_loss作为响应变量，使用性别和锻炼作为两个预测变量。

然后我们可以使用summary()函数来显示模型的结果：

 #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

从模型结果中我们可以看到，性别、运动以及两个变量之间的交互作用在 0.05 显着性水平上均具有统计显着性。

检查模型假设

在进一步进行之前，我们需要验证模型的假设是否得到满足，以便模型结果可靠。特别是，双向方差分析假设：

1. 独立性——每组的观察结果必须相互独立。由于我们采用的是随机设计，所以这个假设应该是满足的，所以我们不必太担心。

2. 正态性——对于两个因素组的每个组合，因变量应具有近似正态分布。

测试此假设的一种方法是创建模型残差的直方图。如果残差近似正态分布，则应满足该假设。

 #define model residuals
reside <- model$residuals

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

残差近似正态分布，因此我们可以假设满足正态性假设。

3. 方差相等——每组的方差相等或近似相等。

检查此假设的一种方法是使用car包执行 Levene 方差齐性检验：

 #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  

由于检验的 p 值大于我们的显着性水平 0.05，因此我们可以假设满足组间方差相等的假设。

分析治疗差异

一旦我们验证了模型假设得到满足，我们就可以进行事后测试来准确确定哪些治疗组彼此不同。

对于我们的事后测试，我们将使用TukeyHSD()函数执行 Tukey 测试以进行多重比较：

 #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

p 值表示每组之间是否存在统计显着差异。

例如，在上面的最后一行中，我们看到，与不运动的女性组相比，不运动的男性组在体重减轻方面没有出现统计学上的显着差异（p 值：0.990364）。

我们还可以使用 R 中的plot()函数可视化 Tukey 测试产生的 95% 置信区间：

 #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)

报告双向方差分析结果

最后，我们可以用总结结果的方式报告双向方差分析的结果：

进行双向方差分析来检查性别（男性、女性）和运动计划（无、轻度、剧烈）对体重减轻（以磅为单位）的影响。性别和运动对减肥的影响之间存在统计学上显着的交互作用（F(2, 54) = 4.615，p = 0.0141）。进行事后 Tukey 的 HSD 测试。

对于男性来说，剧烈运动计划比轻度运动计划 (p < 0.0001) 或无运动计划(p < 0.0001) 能显着减轻体重。此外，对于男性来说，清淡饮食比不运动的减肥效果显着更大（p < 0.0001）。

对于女性来说，剧烈运动计划比轻度运动计划 (p < 0.0001) 或无运动计划(p < 0.0001) 带来的减肥效果显着更大。

进行正态性检查和 Levene 检验以验证方差分析的假设是否得到满足。