Как провести ANCOVA в R

В этом руководстве представлен пример выполнения ANCOVA в R.

Пример: ANCOVA в R

Мы проведем 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)

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

# technique 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(technique) %>%
 summarise (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

Мы видим, что среднее значение и стандартное отклонение текущей оценки для студентов, использующих каждый метод обучения, примерно одинаковы.

Мы также можем видеть, что средний экзаменационный балл заметно выше у студентов, которые использовали технику обучения А , по сравнению с методиками В и С.

Мы также можем визуализировать распределение экзаменационных баллов в зависимости от техники обучения с помощью диаграмм :

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

Проверка предположений ANCOVA с помощью коробчатых диаграмм

Точно так же мы можем использовать диаграммы для визуализации распределения текущих оценок в зависимости от техники обучения:

boxplot(current_grade ~ technique,
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)
#technique 2 74 37.21 0.599 0.552
#Residuals 87 5406 62.14

Значение p больше 0,05, поэтому ковариата ( текущая оценка) и лечение ( техника изучения ) кажутся независимыми.

Затем, чтобы убедиться, что среди групп существует однородность дисперсии, мы можем провести тест Левена:

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

#conduct Levene's Test
leveneTest(exam~technique, 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() в пакете car, чтобы указать, что мы хотели бы использовать для модели сумму квадратов типа 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") 

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

Мы видим, что значение p для техники чрезвычайно мало, что указывает на то, что изучение техники оказывает статистически значимое влияние на результаты экзамена, даже с поправкой на текущую оценку.

Шаг 4: Апостериорные тесты

Хотя результаты ANCOVA показали нам, что техника обучения оказывает статистически значимое влияние на результаты экзаменов, нам необходимо провести апостериорные тесты , чтобы выяснить, какие методы обучения отличаются друг от друга.

Для этого мы можем использовать функцию glht() в пакете multcomp в R для выполнения теста Тьюки для множественных сравнений:

#load the *multcomp* library
library(multcomp)

#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 Hypotheses:
# 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 Hypotheses:
# 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

Из вывода мы видим, что существует статистически значимая разница (при α = 0,05) в экзаменационных баллах между изучением техники A и техникой B (значение p: 0,000109), а также между техникой A и техникой C ( p-значение: <1e-04).

Мы также можем видеть, что нет статистически значимой разницы (при α = 0,05) между методами B и C. Доверительные интервалы между методиками также подтверждают эти выводы.

Таким образом, можно сделать вывод, что использование методики обучения А приводит к статистически значимому увеличению экзаменационного балла учащихся по сравнению с методиками В и С , даже с поправкой на текущую оценку учащегося в классе.