如何在 r 中进行 ancova

本教程提供了如何在 R 中执行ANCOVA的示例。

示例：R 中的 ANCOVA

我们将使用以下变量执行 ANCOVA 来测试学习技术是否对考试结果产生影响：

• 技术研究：我们希望分析的自变量
• 当前学生成绩：我们希望考虑的协变量
• 评论分数：我们要分析的响应变量

以下数据集包含 90 名学生的信息，随机分为三组，每组 30 人。

该数据集显示了每个学生使用的学习技巧（A、B或C） 、他们开始使用该学习技巧时在班级中的当前成绩，以及使用该学习技巧准备一个月后在考试中获得的成绩为了考试。考试：

 #make this example reproducible
set.seed(10)

#create dataset
data <- data.frame(technique = rep(c("A", "B", "C"), each = 30),
                   current_grade = runif(90, 65, 95),
                   exam = c(runif(30, 80, 95), runif(30, 70, 95), runif(30, 70, 90)))

#view first six lines of dataset
head(data)

# technical current_grade exam
#1 A 80.22435 87.32759
#2 A 74.20306 90.67114
#3 A 77.80723 88.87902
#4 A 85.79306 87.75735
#5 A 67.55408 85.72442
#6 A 71.76310 92.52167

第 1 步：探索数据

在拟合 ANCOVA 模型之前，我们首先需要探索数据以更好地理解它并验证不存在可能扭曲结果的极端异常值。

首先，我们可以显示数据集中每个变量的摘要：

 summary(data)

# technical current_grade exam      
#A:30 Min. :65.43 Min. :71.17  
# B:30 1st Qu.:71.79 1st Qu.:77.27  
# C:30 Median:77.84 Median:84.69  
# Mean:78.15 Mean:83.38  
# 3rd Qu.:83.65 3rd Qu.:89.22  
# Max. :93.84 Max. :94.76  

我们可以看到每个研究技术值（ A、B和C）在数据中出现了 30 次。

我们还可以看到学生当前分数在研究开始时的分布情况。全班最低分65.43分，最高分93.84分，平均分78.15分。

同样，我们可以看到，这次考试获得的最低分数是71.17，最高分数是94.76，平均分数是83.38。

然后我们可以使用dplyr包轻松找到每种学习技术当前成绩和考试结果的平均值和标准差：

 #load dplyr
library(dplyr)

data %>%
  group_by (technical) %>%
  summarize (mean_grade = mean(current_grade),
            sd_grade = sd(current_grade),
            mean_exam = mean(exam),
            sd_exam = sd(exam))

# A tibble: 3 x 5
# technique mean_grade sd_grade mean_exam sd_exam                      
#1 A 79.0 7.00 88.5 3.88
#2 B 78.5 8.33 81.8 7.62
#3 C 76.9 8.24 79.9 5.71

我们可以看到，使用每种学习技术的学生当前年级的平均值和标准差大致相似。

我们还可以看到，与技术B和C相比，使用学习技术A的学生的平均考试成绩明显更高。

我们还可以使用箱线图根据学习技术可视化考试结果的分布：

 boxplot(exam ~ technique,
data = data,
main = "Exam Score by Studying Technique",
xlab = "Studying Technique",
ylab = "Exam Score",
col = "steelblue",
border = "black"
)

同样，我们还可以使用箱线图根据学习技术来可视化当前成绩的分布：

 boxplot(current_grade ~ technical,
data = data,
main = "Current Grade by Studying Technique",
xlab = "Studying Technique",
ylab = "Current Grade",
col = "steelblue",
border = "black"
)

第 2 步：检查模型假设

一旦我们进行了一些基本的数据探索并熟悉了数据，我们就需要验证是否满足以下 ANCOVA 假设：

• 协变量和治疗是独立的——有必要验证协变量（当前成绩）和治疗（研究技术）是相互独立的，因为只有当协变量和治疗独立于反应变量（检查）。
• 方差同质性——我们需要检查组之间的方差是否相等

为了验证协变量和治疗是独立的，我们可以使用当前成绩作为响应变量并使用研究技术作为预测变量来执行方差分析：

 #fit anova model
anova_model <- aov(current_grade ~ technique, data = data)
#view summary of anova model
summary(anova_model)

# Df Sum Sq Mean Sq F value Pr(>F)
#technical 2 74 37.21 0.599 0.552
#Residuals 87 5406 62.14    

p 值大于 0.05，因此协变量（当前等级）和治疗（研究技术）显得独立。

然后，为了检查组之间是否存在方差同质性，我们可以执行 Levene 检验：

 #load car library to conduct Levene's Test
libary(car)

#conduct Levene's Test
leveneTest(technical exam, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)    
#group 2 9.4324 0.0001961 ***
#87   

检验的 p 值等于 0.0001961，这表明组之间的方差不相等。尽管我们可以尝试对数据进行转换来纠正这个问题，但此时我们不会太担心方差差异。

步骤 3：调整 ANCOVA 模型

接下来，我们将使用考试成绩作为响应变量、学习技术作为预测（或“治疗”）变量以及当前成绩作为协变量来拟合 ANCOVA 模型。

我们将使用包中的 Anova() 函数，因为要做到这一点，只是为了能够指定我们想要对模型使用 III 型平方和，因为 I 型平方和取决于将预测变量输入到模型中：

 #load car library
library(car)

#fit ANCOVA model
ancova_model <- aov(exam ~ technique + current_grade, data = data)

#view summary of model
Anova(ancova_model, type="III") 

#Answer: exam
# Sum Sq Df F value Pr(>F)    
#(Intercept) 7161.2 1 201.4621 < 2.2e-16 ***
#technical 1242.9 2 17.4830 4.255e-07 ***
#current_grade 12.3 1 0.3467 0.5576    
#Residuals 3057.0 86         

我们可以看到，技术的 p 值极低，这表明即使在控制当前成绩之后，学习技术对考试成绩也有统计上的显着影响。

第 4 步：事后测试

尽管 ANCOVA 结果告诉我们，学习技术对考试成绩有统计上的显着影响，但我们需要进行事后测试，以确定哪些学习技术彼此不同。

为此，我们可以使用 R 中的multcomp包中的 glht() 函数来执行 Tukey 测试以进行多重比较：

 #load the multcomp library
library(multicomp)

#fit the ANCOVA model
ancova_model <- aov(exam ~ technique + current_grade, data = data)

#define the post hoc comparisons to make
postHocs <- glht(ancova_model, linfct = mcp(technique = "Tukey"))

#view a summary of the post hoc comparisons
summary(postHocs)

#Multiple Comparisons of Means: Tukey Contrasts
#
#Fit: aov(formula = exam ~ technique + current_grade, data = data)
#
#Linear Assumptions:
#Estimate Std. Error t value Pr(>|t|)    
#B - A == 0 -6.711 1.540 -4.358 0.000109 ***
#C - A == 0 -8.736 1.549 -5.640 < 1e-04 ***
#C - B == 0 -2.025 1.545 -1.311 0.393089    

#view the confidence intervals associated with the multiple comparisons
confint(postHocs)

# Simultaneous Confidence Intervals
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#Fit: aov(formula = exam ~ technique + current_grade, data = data)
#
#Quantile = 2.3845
#95% family-wise confidence level
#
#Linear Assumptions:
# Estimate lwr upr     
#B - A == 0 -6.7112 -10.3832 -3.0392
#C - A == 0 -8.7364 -12.4302 -5.0426
#C - B == 0 -2.0252 -5.7091 1.6588

从结果中，我们可以看到，学习技术A和学习技术B （p 值：0,000109）以及技术A和技术C之间的考试成绩存在统计显着性差异（α = 0.05） （p 值：<1e-04）。

我们还可以看到，技术B和C之间没有统计学上的显着差异（α = 0.05）。技术之间的置信区间也支持这些发现。

因此，我们可以得出结论，即使在控制了学生当前在班级中的成绩之后，与技术B和C相比，使用学习技术A会导致学生的考试成绩在统计上显着提高。