# 安裝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語言和統計首發，轉載請註明
▌編輯：June