Як провести 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)

# 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

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

Ми також бачимо, що середній бал на іспиті значно вищий для студентів, які використовували метод навчання A , порівняно з методами B і C.

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

 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, тому коваріата ( поточна оцінка) і лікування ( техніка дослідження ) здаються незалежними.

Потім, щоб перевірити наявність однорідності дисперсії між групами, ми можемо виконати тест Левена:

 #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 показали нам, що методика дослідження мала статистично значущий вплив на оцінки іспитів, нам потрібно провести нестандартне тестування , щоб визначити, які методи дослідження відрізняються одна від одної.

Для цього ми можемо використати функцію glht() у пакеті multcomp у R, щоб виконати тест 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

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

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

Таким чином, ми можемо зробити висновок, що використання методу навчання A призводить до статистично значуще вищої оцінки іспиту для студентів порівняно з методами B і C , навіть після контролю поточної оцінки студента в класі.