R語言統計篇: 單因素協方差分析

單因素協方差分析（One-way ANCOVA）可以研究一個分類變量對一個連續變量的影響，同時校正其他變量的作用，這些變量也稱為協變量（Covariate）。也是單因素方差分析（One-way ANOVA，R語言統計篇：單因素方差分析）的一個延伸。比方說，我們現在想要研究不同BMI（偏輕，正常與超重）與空腹血糖的關係，同時校正血壓水平。在此研究中，BMI分組是一個分類變量（自變量），血糖是一個連續變量（因變量），血壓則是一個協變量（covariate）。c. 在各個組別中，因變量與協變量之間需呈線性關係；d. 在各個組別中，因變量與協變量所形成的回歸直線的斜率相等。

協方差分析的理論知識較複雜，感興趣的讀者可以參考其他教材讀物。本公眾號更偏重實戰，直接分享代碼迅速帶你走起來。使用到的數據來自{datarium}裡面的「stress」，首先安裝{dataruim}。

# 安裝datarium
install.packages("datarium")
# 載入
library(datarium)
# 查看stress的基本信息
summary(stress)
 #      id            score        treatment     exercise       age       
 #Min.   : 1.00   Min.   : 65.80   yes:30    low     :20   Min.   :52.00  
 #1st Qu.:15.75   1st Qu.: 80.40   no :30    moderate:20   1st Qu.:57.00  
 #Median :30.50   Median : 85.55             high    :20   Median :60.00  
 #Mean   :30.50   Mean   : 84.58                           Mean   :59.95  
 #3rd Qu.:45.25   3rd Qu.: 90.80                           3rd Qu.:62.00  
 #Max.   :60.00   Max.   :100.00                           Max.   :75.00
「stress」數據集記錄了幾個變量與壓力水平（stress）的關係。為了簡化上述的數據，小編只提取三個變量：score，exercise與age。於是，我們的研究問題變成了：校正協變量年齡（age）後，運動（exercise）是否與壓力水平（score）有關？下一步，提取三個變量（age, exercise與score）： 
# 提取三個變量
mydata <- stress[, c("score", "exercise", "age")]
# 查看mydata
summary(mydata)
 #    score            exercise       age       
 #Min.   : 65.80   low     :20   Min.   :52.00  
 #1st Qu.: 80.40   moderate:20   1st Qu.:57.00  
 #Median : 85.55   high    :20   Median :60.00  
 #Mean   : 84.58                 Mean   :59.95  
 #3rd Qu.: 90.80                 3rd Qu.:62.00  
 #Max.   :100.00                 Max.   :75.00

其中，exercise包含三組：low, moderate與high，分別代表運動強度的低、中、高組。
# 載入car，需提前安裝
library(car)
# Levene檢驗
leveneTest(mydata$score, mydata$exercise)
#Levene's Test for Homogeneity of Variance (center = median)
#      Df F value Pr(>F)
#group  2  0.9593 0.3893
#      57
其中，Pr(>F) 為0.3893，提示各組別間的方差沒有統計學差異，即未違反方差齊性原則。
# 建立回歸模型，協變量age需在exercise前
myfit1 <- lm(score ~ age + exercise, data = mydata)
# 將殘差作直方圖
hist(residuals(myfit1))
2.3 在各個組別中，因變量與協變量之間需呈線性關係
library(ggplot2)
ggplot(mydata, aes(age, score, color = exercise)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = F)
通過肉眼觀察，age與score在各組間大體上呈線性關係。在回歸方程中添加交互項，如果交互項存在顯著性，那麼就違反此原則；反之，則沒有違反。
# 含有交互項
myfit2 <- lm(score ~ age + exercise + age:exercise, data = mydata)
Anova(myfit2, type = "III")
#Anova Table (Type III tests)

#Response: score
#              Sum Sq Df F value   Pr(>F)   
#(Intercept)   374.13  1 11.2246 0.001478 **
#age            65.62  1  1.9686 0.166319   
#exercise       21.05  2  0.3158 0.730555   
#age:exercise   13.70  2  0.2055 0.814867   
#Residuals    1799.88 54                    
#---
#Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

myfit2中的age:exercise即為交互項，從結果可知，交互項的p值為0.815，提示沒有統計學意義，即沒有違反此前提。由於我們沒有違反協方差分析的前提，接下來就可以進行協方差分析了，並且毫無後顧之憂。
# 協方差分析
myfit1 <- lm(score ~ age + exercise, data = mydata)
Anova(myfit1, type = "III")
#Anova Table (Type III tests)

#Response: score
#             Sum Sq Df F value    Pr(>F)    
#(Intercept)  670.77  1 20.7121 2.917e-05 ***
#age          298.42  1  9.2147  0.003639 ** 
#exercise     980.91  2 15.1444 5.528e-06 ***
#Residuals   1813.58 56                      
#---
#Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
從上述結果中可以得出以下結論：在校正年齡後，運動與壓力水平有關，並且結果存在統計學意義（p = 5.528e-06）。
目前只知道運動與壓力水平有關，但不知道具體哪些組之間存在差異，因此需要進行多重比較來回答這個問題。
# 安裝所需R包
install.packages("emmeans")
install.packages("rstatix")

# 載入
library(emmeans)
library(rstatix)

# 兩兩比較
pwc <- emmeans_test(
    score ~ exercise, covariate = age,
    p.adjust.method = "bonferroni", 
    data = mydata
    )
pwc
## A tibble: 3 x 9
#  term         .y.   group1   group2      df statistic          p     p.adj p.adj.signif
#* <chr>        <chr> <chr>    <chr>    <dbl>     <dbl>      <dbl>     <dbl> <chr>       
#1 age*exercise score low      moderate    56    -0.105 0.917      1         ns          
#2 age*exercise score low      high        56     4.68  0.0000186  0.0000557 ****        
#3 age*exercise score moderate high        56     5.02  0.00000558 0.0000167 ****

# A tibble: 3 x 8
    age exercise emmean    se    df conf.low conf.high method      
  <dbl> <fct>     <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
1  60.0 low        87.6  1.32    56     85.0      90.3 Emmeans test
2  60.0 moderate   87.8  1.28    56     85.2      90.4 Emmeans test
3  60.0 high       78.3  1.36    56     75.6      81.0 Emmeans test
high組（78.3±1.36）的壓力水平要顯著低於low組（87.6± 1.32）與moderate組（87.8 ± 1.28），而low組與moderate組沒有統計學差異。

Discovering statistics using R, Andy Field.

讓R語言與統計變得簡單！