Cách thực hiện kiểm định ANCOVA trong R
Duc Nguyen | Chuyên đào tạo kỹ năng R |
www.tuhocr.com
August 3, 2023
Phân biệt giữa ANCOVA và ANOVA
Phân tích hiệp phương sai Analysis of Covariance (ANCOVA) là trường hợp mở rộng của phân tích phương sai Analysis of Variance (ANOVA) khi cho phép phân tích thêm biến C (continuous) tác động lên biến Y (continuous) để tìm ra sự khác biệt giữa các nhóm trong biến X (category).
Ứng dụng dễ thấy nhất là nếu thực hiện phân tích phương sai ANOVA 1 yếu tố thông thường, giữa các nhóm trong biến X tác động lên biến Y mà có p-value > 0.05 chẳng hạn, tức là không thỏa điều kiện về ý nghĩa thống kê, thì ta có thể đánh giá thêm xem biến Y bản chất có bị ảnh hưởng bởi 1 biến C nào đó (ở dạng continuous) hay không, khi đó ta đưa biến C vào để hiệu chỉnh lại kết quả từ Y, rồi thực hiện phép phân tích ANCOVA để tìm ra p-value lúc này ở dạng kết quả hiệu chỉnh cho phép đánh giá chính xác hơn sự khác biệt giữa các nhóm trong biến X (category).
ANCOVA có hai loại one-way và two-way tương tự như ANOVA, tương ứng cho số lượng yếu tố X trong mô hình.
Các giả định của ANCOVA
Linearity between the covariate and the outcome variable at each level of the grouping variable. This can be checked by creating a grouped scatter plot of the covariate and the outcome variable.
Homogeneity of regression slopes. The slopes of the regression lines, formed by the covariate and the outcome variable, should be the same for each group. This assumption evaluates that there is no interaction between the outcome and the covariate. The plotted regression lines by groups should be parallel.
The outcome variable should be approximately normally distributed. This can be checked using the Shapiro-Wilk test of normality on the model residuals.
Homoscedasticity or homogeneity of residuals variance for all groups. The residuals are assumed to have a constant variance (homoscedasticity)
No significant outliers in the groups
Case study cho one-way ANCOVA
Researchers investigated the effect of exercises in reducing the level of anxiety. Therefore, they conducted an experiment, where they measured the anxiety score of three groups of individuals practicing physical exercises at different levels (grp1: low, grp2: moderate and grp3: high). The anxiety score was measured pre- and 6-months post-exercise training programs. It is expected that any reduction in the anxiety by the exercises programs would also depend on the participant’s basal level of anxiety score. In this analysis we use the pretest anxiety score as the covariate and are interested in possible differences between group with respect to the post-test anxiety scores.
Tạm dịch: Nghiên cứu này từ dataset anxiety của package
datarium. Mục tiêu của nghiên cứu là khảo sát ảnh hưởng của
việc tập thể dục lên khả năng giảm căng thẳng của sinh viên.
Trong đó theo dõi 45 sinh viên, chia làm 3 nhóm (15 bạn/nhóm) ở 3 mức
độ khác nhau về việc tập thể dục, ký hiệu tương ứng là grp1
tập thể dục ít, grp2 tập thể dục trung bình,
grp3 tập thể dục nhiều. Sau đó các nhà nghiên cứu đo điểm
căng thẳng của từng bạn sinh viên xem như thế nào ở 6 tháng trước khi
bắt đầu tập thể dục (pretest) và 6 tháng sau khi tập thể
dụng (postest).
Câu hỏi đặt ra là: kết quả của việc tập thể dục theo các cường độ khác nhau như vậy (ít, vừa, nhiều) có làm giảm căng thẳng hay không?
Nếu bình thường chỉ phân tích ANOVA thì ta chỉ lấy biến
postest là biến y để làm cơ sở đưa vào phân tích one-way
ANOVA với yếu tố là mức độ tập thể dục. Tuy nhiên rõ ràng giữa các bạn
sinh viên này thì trước khi tham gia chương trình tập thể
dục này các bạn đã có mức độ căng thẳng khác nhau rồi (đo thông qua
pretest) nên nếu chúng ta không hiệu chỉnh
giá trị này thì kết quả thu được từ phân tích phương sai sẽ bị ảnh hưởng
(lỡ nhóm tập thể dục nhiều, theo lý thuyết sẽ giảm căng thẳng, nhưng vì
các bạn sinh viên chọn ngẫu nhiên trong nhóm này có sự căng thẳng
nền trước đó dẫn đến cho dù có tập thể dục
nhiều thì điểm căng thẳng cũng không giảm). Đây là lý do vì sao chúng ta
cần đưa covariate (gọi là hiệp biến) ở đây là yếu tố căng thẳng nền vào,
hoặc có thể đưa thêm yếu tố độ tuổi vào cũng được, nghĩa là covariate có
thể từ 1 trở lên để tác động đến biến y sau cùng, giúp điều chỉnh kết
quả của y (postest) tương ứng với từng nhóm trong biến x (mức độ tập thể
dục) được chính xác hơn.
Dữ liệu cụ thể
anxiety <- read.csv("ancova_example_1.csv")
anxiety ## id group pretest posttest
## 1 1 grp1 14.1 14.1
## 2 2 grp1 14.5 14.3
## 3 3 grp1 15.7 14.9
## 4 4 grp1 16.0 15.3
## 5 5 grp1 16.5 15.7
## 6 6 grp1 16.9 16.2
## 7 7 grp1 17.0 16.5
## 8 8 grp1 17.0 16.6
## 9 9 grp1 17.3 16.5
## 10 10 grp1 17.3 16.7
## 11 11 grp1 17.8 17.3
## 12 12 grp1 17.9 17.5
## 13 13 grp1 19.1 19.3
## 14 14 grp1 19.4 19.0
## 15 15 grp1 19.8 19.4
## 16 16 grp2 13.7 12.7
## 17 17 grp2 14.7 13.1
## 18 18 grp2 14.9 13.6
## 19 19 grp2 15.1 13.6
## 20 20 grp2 15.8 14.2
## 21 21 grp2 16.4 14.9
## 22 22 grp2 16.6 16.1
## 23 23 grp2 16.9 16.1
## 24 24 grp2 16.9 16.3
## 25 25 grp2 17.2 15.9
## 26 26 grp2 17.8 17.4
## 27 27 grp2 17.8 16.9
## 28 28 grp2 18.2 17.1
## 29 29 grp2 18.4 17.3
## 30 30 grp2 19.3 17.7
## 31 31 grp3 14.6 11.7
## 32 32 grp3 15.0 11.9
## 33 33 grp3 15.5 11.0
## 34 34 grp3 15.7 12.1
## 35 35 grp3 16.4 12.3
## 36 36 grp3 16.9 13.6
## 37 37 grp3 17.1 14.3
## 38 38 grp3 17.3 14.2
## 39 39 grp3 17.5 14.4
## 40 40 grp3 17.6 13.8
## 41 41 grp3 17.8 14.3
## 42 42 grp3 17.9 13.8
## 43 43 grp3 18.4 15.4
## 44 44 grp3 18.5 15.1
## 45 45 grp3 19.0 15.5Kiểm tra giả thuyết thống kê
1. Linearity assumption
There was a linear relationship between pre-test and post-test anxiety score for each training group, as assessed by visual inspection of a scatter plot.
library(tidyverse)
library(ggpubr)
library(rstatix)
library(broom)
ggscatter(data = anxiety, x = "pretest", y = "posttest", color = "group", add = "reg.line") +
stat_regline_equation(
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~~~~~~"),
color = group),
show.legend = NULL
)
2. Homogeneity of regression slopes
This assumption checks that there is no significant interaction between the covariate and the grouping variable.
There was homogeneity of regression slopes as the interaction term was not statistically significant, F(2, 39) = 0.13, p = 0.88.
anxiety %>% anova_test(posttest ~ group*pretest)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 group 2 39 209.314 1.40e-21 * 0.915
## 2 pretest 1 39 572.828 6.36e-25 * 0.936
## 3 group:pretest 2 39 0.127 8.81e-01 0.0063. Normality of residuals
# Fit the model, the covariate goes first
model <- lm(posttest ~ pretest + group, data = anxiety)
# Inspect the model diagnostic metrics
model.metrics <- augment(model)[, c(1, 2, 3, 5, 8, 9)]
model.metrics <- as.data.frame(model.metrics)
model.metrics## posttest pretest group .resid .cooksd .std.resid
## 1 14.1 14.1 grp1 0.549556955 1.008643e-01 1.457068239
## 2 14.3 14.5 grp1 0.338455577 3.095139e-02 0.885202791
## 3 14.9 15.7 grp1 -0.294848556 1.334801e-02 -0.749866153
## 4 15.3 16.0 grp1 -0.203174590 5.675944e-03 -0.514543486
## 5 15.7 16.5 grp1 -0.317051312 1.206505e-02 -0.799162006
## 6 16.2 16.9 grp1 -0.228152690 5.919695e-03 -0.574088378
## 7 16.5 17.0 grp1 -0.030928035 1.082522e-04 -0.077810670
## 8 16.6 17.0 grp1 0.069071965 5.399278e-04 0.173775537
## 9 16.5 17.3 grp1 -0.339254068 1.311360e-02 -0.853697710
## 10 16.7 17.3 grp1 -0.139254068 2.209467e-03 -0.350418434
## 11 17.3 17.8 grp1 -0.053130791 3.483819e-04 -0.134045731
## 12 17.5 17.9 grp1 0.044093865 2.461007e-04 0.111341667
## 13 19.3 19.1 grp1 0.610789731 7.558967e-02 1.572623744
## 14 19.0 19.4 grp1 0.002463698 1.425655e-06 0.006392114
## 15 19.4 19.8 grp1 -0.008637680 2.154133e-05 -0.022683279
## 16 12.7 13.7 grp2 0.201780151 1.330855e-02 0.534203822
## 17 13.1 14.7 grp2 -0.425973294 3.561620e-02 -1.095079357
## 18 13.6 14.9 grp2 -0.131523983 3.095512e-03 -0.336670267
## 19 13.6 15.1 grp2 -0.337074672 1.862872e-02 -0.859563007
## 20 14.2 15.8 grp2 -0.456502083 2.661806e-02 -1.153075147
## 21 14.9 16.4 grp2 -0.373154150 1.590854e-02 -0.939084242
## 22 16.1 16.6 grp2 0.621295161 4.364296e-02 1.563046354
## 23 16.1 16.9 grp2 0.312969127 1.119730e-02 0.787636564
## 24 16.3 16.9 grp2 0.512969127 3.008104e-02 1.290968359
## 25 15.9 17.2 grp2 -0.195356906 4.550936e-03 -0.492313232
## 26 17.4 17.8 grp2 0.687991027 6.659086e-02 1.743825167
## 27 16.9 17.8 grp2 0.187991027 4.971911e-03 0.476493837
## 28 17.1 18.2 grp2 -0.023110351 8.781534e-05 -0.058939936
## 29 17.3 18.4 grp2 -0.028661040 1.474391e-04 -0.073375382
## 30 17.7 19.3 grp2 -0.553639140 8.575232e-02 -1.451061685
## 31 11.7 14.6 grp3 0.620311647 9.507982e-02 1.613950448
## 32 11.9 15.0 grp3 0.409210269 3.392908e-02 1.053609373
## 33 11.0 15.5 grp3 -1.004666453 1.631872e-01 -2.560468115
## 34 12.1 15.7 grp3 -0.110217142 1.812547e-03 -0.279990633
## 35 12.3 16.4 grp3 -0.629644554 4.784434e-02 -1.587371864
## 36 13.6 16.9 grp3 0.156478724 2.773679e-03 0.393690617
## 37 14.3 17.1 grp3 0.650928035 4.795104e-02 1.637645152
## 38 14.2 17.3 grp3 0.345377346 1.367987e-02 0.869284480
## 39 14.4 17.5 grp3 0.339826657 1.360708e-02 0.856054572
## 40 13.8 17.6 grp3 -0.362948688 1.581104e-02 -0.914851288
## 41 14.3 17.8 grp3 -0.068499377 5.897635e-04 -0.172926642
## 42 13.8 17.9 grp3 -0.671274721 5.820829e-02 -1.696230176
## 43 15.4 18.4 grp3 0.414848556 2.642395e-02 1.055053126
## 44 15.1 18.5 grp3 0.012073212 2.330621e-05 0.030755387
## 45 15.5 19.0 grp3 -0.101803510 2.073307e-03 -0.261953948# Assess normality of residuals using shapiro wilk test
shapiro_test(model.metrics$.resid)## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 model.metrics$.resid 0.975 0.444The Shapiro Wilk test was not significant (p > 0.05), so we can assume normality of residuals.
4. Homogeneity of variances
ANCOVA assumes that the variance of the residuals is equal for all groups. This can be checked using the Levene’s test.
The Levene’s test was not significant (p > 0.05), so we can assume homogeneity of the residual variances for all groups.
model.metrics %>% levene_test(.resid ~ group)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 42 2.27 0.1165. Outliers
There were no outliers in the data, as assessed by no cases with standardized residuals greater than 3 in absolute value.
model.metrics %>%
filter(abs(.std.resid) > 3) %>%
as.data.frame()## [1] posttest pretest group .resid .cooksd .std.resid
## <0 rows> (or 0-length row.names)Tính toán ANCOVA
Sau khi dataset đã thỏa các điều kiện để thực hiện ANCOVA, giờ ta đưa vào model để tính ra kết quả.
After adjustment for pre-test anxiety score, there was a statistically significant difference in post-test anxiety score between the groups, F(2, 41) = 218.63, p < 0.0001.
res.aov <- anxiety %>% anova_test(posttest ~ pretest + group) # pretest là covariate phải đặt phía trước.
get_anova_table(res.aov)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 pretest 1 41 598.321 4.48e-26 * 0.936
## 2 group 2 41 218.629 1.35e-22 * 0.914Thực hiện phân hạng Post-hoc test
Pairwise comparisons can be performed to identify which groups are different. The Bonferroni multiple testing correction is applied.
# Pairwise comparisons
library(emmeans)
pwc <- anxiety %>%
emmeans_test(
posttest ~ group, covariate = pretest,
p.adjust.method = "bonferroni"
)
pwc## # A tibble: 3 × 9
## term .y. group1 group2 df statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 pretest*gr… post… grp1 grp2 41 4.24 1.26e- 4 3.77e- 4 ***
## 2 pretest*gr… post… grp1 grp3 41 19.9 1.19e-22 3.58e-22 ****
## 3 pretest*gr… post… grp2 grp3 41 15.5 9.21e-19 2.76e-18 ****# Display the adjusted means of each group
# Also called as the estimated marginal means (emmeans)
get_emmeans(pwc)## # A tibble: 3 × 8
## pretest group emmean se df conf.low conf.high method
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 16.9 grp1 16.4 0.106 41 16.2 16.7 Emmeans test
## 2 16.9 grp2 15.8 0.107 41 15.6 16.0 Emmeans test
## 3 16.9 grp3 13.5 0.106 41 13.2 13.7 Emmeans testTrình bày kết quả
An ANCOVA was run to determine the effect of exercises on the anxiety score after controlling for basal anxiety score of participants.
After adjustment for pre-test anxiety score, there was a statistically significant difference in post-test anxiety score between the groups, F(2, 41) = 218.63, p < 0.0001.
Post hoc analysis was performed with a Bonferroni adjustment. The mean anxiety score was statistically significantly greater in grp1 (16.4 +/- 0.15) compared to the grp2 (15.8 +/- 0.12) and grp3 (13.5 +/_ 0.11), p < 0.001.
Tạm dịch: Sau khi thực hiện phép phân tích ANCOVA với
việc hiệu chỉnh biến pretest thì ta thấy kết quả tập thể
dục của các nhóm thật sự khác biệt nhau, trong đó nhóm tập thể dục nhiều
(grp3) có chỉ số căng thẳng thấp nhất.
# Visualization: line plots with p-values
pwc <- pwc %>% add_xy_position(x = "group", fun = "mean_se")
ggline(get_emmeans(pwc), x = "group", y = "emmean") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)