วิธีดำเนินการ ancova ใน r

บทช่วยสอนนี้ให้ตัวอย่างวิธีดำเนินการ ANCOVA ใน R

ตัวอย่าง: ANCOVA ใน R

เราจะทำ ANCOVA เพื่อทดสอบว่าเทคนิคการเรียนมีผลกระทบต่อผลการสอบหรือไม่โดยใช้ตัวแปรต่อไปนี้:

• การศึกษาทางเทคนิค : ตัวแปรอิสระที่เราต้องการวิเคราะห์
• เกรดของนักเรียนปัจจุบัน : ตัวแปรร่วมที่เราต้องการนำมาพิจารณา
• คะแนนรีวิว : ตัวแปรการตอบสนองที่เราต้องการวิเคราะห์

ชุดข้อมูลต่อไปนี้ประกอบด้วยข้อมูลเกี่ยวกับนักเรียน 90 คน สุ่มแบ่งออกเป็น 3 กลุ่ม กลุ่มละ 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

นอกจากนี้เรายังสามารถแสดงภาพการกระจายผลการสอบตามเทคนิคการเรียนโดยใช้ boxplots :

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

ในทำนองเดียวกัน เรายังสามารถใช้ boxplots เพื่อแสดงภาพการแจกแจง เกรดปัจจุบัน ตามเทคนิคการศึกษา:

 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 จะบอกเราว่า เทคนิคการศึกษา มีผลกระทบที่มีนัยสำคัญทางสถิติต่อคะแนนการสอบ แต่เราจำเป็นต้องดำเนิน การทดสอบภายหลัง เพื่อพิจารณาว่าเทคนิคการศึกษาใดที่แตกต่างกัน

เมื่อต้องการทำเช่นนี้ เราสามารถใช้ฟังก์ชัน 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 แม้ว่าหลังจากควบคุมเกรดปัจจุบันของนักเรียนในชั้นเรียนแล้วก็ตาม