前面介绍了多个样本均数的多重比较，多样本非参数检验后的多重比较：
R语言多个样本均数的多重比较
R语言非参数检验后的多重比较
今天学习下重复测量数据的多重比较，本篇内容和课本结果差异较大，如有错误欢迎指出。
使用的数据来自孙振球，徐勇勇《医学统计学》第4版。课本的电子版已上传到QQ群，加群即可免费获取！
课本封面
重复测量方差分析
使用课本例12-1的数据，直接读取：
代码语言：javascript复制df12_3 <- foreign::read.spss("E:/各科资料/医学统计学/研究生课程/析因设计重复测量/9重复测量18-9研/例12-03.sav",to.data.frame = T)
str(df12_3)
## 'data.frame': 15 obs. of 7 variables:
## $ No : num 1 2 3 4 5 6 7 8 9 10 ...
## $ group: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...
## $ t0 : num 120 118 119 121 127 121 122 128 117 118 ...
## $ t1 : num 108 109 112 112 121 120 121 129 115 114 ...
## $ t2 : num 112 115 119 119 127 118 119 126 111 116 ...
## $ t3 : num 120 126 124 126 133 131 129 135 123 123 ...
## $ t4 : num 117 123 118 120 126 137 133 142 131 133 ...
## - attr(*, "variable.labels")= Named chr [1:7] "xd0xf2xbaxc5" "xd7xe9xb1xf0" "" "" ...
## ..- attr(*, "names")= chr [1:7] "No" "group" "t0" "t1" ...
数据一共7列，第1列是患者编号，第2列是诱导方法（3种），第3-7列是5个时间点的血压。
首先转换数据格式：
代码语言：javascript复制library(reshape2)
df.l <- melt(df12_3, id.vars = c("No","group"),
variable.name = "times",
value.name = "hp")
df.l$No <- factor(df.l$No)
str(df.l)
## 'data.frame': 75 obs. of 4 variables:
## $ No : Factor w/ 15 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ group: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...
## $ times: Factor w/ 5 levels "t0","t1","t2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ hp : num 120 118 119 121 127 121 122 128 117 118 ...
head(df.l)
## No group times hp
## 1 1 A t0 120
## 2 2 A t0 118
## 3 3 A t0 119
## 4 4 A t0 121
## 5 5 A t0 127
## 6 6 B t0 121
进行重复测量方差分析，默认方法不能输出球形检验的结果，所以我更推荐rstatix提供的方法：
# 默认
f <- aov(hp ~ group*times Error(No/times), data = df.l)
summary(f)
##
## Error: No
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 912.2 456.1 5.783 0.0174 *
## Residuals 12 946.5 78.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: No:times
## Df Sum Sq Mean Sq F value Pr(>F)
## times 4 2336.5 584.1 106.6 < 2e-16 ***
## group:times 8 837.6 104.7 19.1 1.62e-12 ***
## Residuals 48 263.1 5.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rstatix
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
anova_test(data = df.l,
dv = hp,
wid = No,
within = times,
between = group
)
## ANOVA Table (type II tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 group 2 12 5.783 1.70e-02 * 0.430
## 2 times 4 48 106.558 3.02e-23 * 0.659
## 3 group:times 8 48 19.101 1.62e-12 * 0.409
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 times 0.293 0.178
## 2 group:times 0.293 0.178
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 times 0.679 2.71, 32.58 1.87e-16 * 0.896 3.59, 43.03 4.65e-21
## 2 group:times 0.679 5.43, 32.58 4.26e-09 * 0.896 7.17, 43.03 2.04e-11
## p[HF]<.05
## 1 *
## 2 *
画图展示：
代码语言：javascript复制library(ggplot2)
df.l |>
group_by(times,group) |>
summarise(mm=mean(hp)) |>
ggplot(aes(times,mm))
geom_line(aes(group=group,color=group),size=1.2)
theme_bw()
接下来是重复测量数据的多重比较，课本中分成了3个方面。
组间差别多重比较
LSD/SNK/Tukey/Dunnett/Bonferroni等方法都可以，和多个均数比较的多重检验一样。
library(PMCMRplus)
summary(lsdTest(hp ~ group, data = df.l))
##
## Pairwise comparisons using Least Significant Difference Test
## data: hp by group
## alternative hypothesis: two.sided
## P value adjustment method: none
## H0
## t value Pr(>|t|)
## B - A == 0 2.175 0.0329218 *
## C - A == 0 3.860 0.0002446 ***
## C - B == 0 1.686 0.0962097 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P值和课本不太一样，但是结论是一样的，A组和B组之间，A组和C组之间有差别，B组和C组之间没有差别。
时间趋势比较
重复测量方差分析可以采取正交多项式来探索时间变化趋势，具体的内涵解读可以参考冯国双老师的这篇文章：https://mp.weixin.qq.com/s/ndinwbDJsHjAelvNfwqgwA
在R里面进行正交多项式的探索略显复杂，首先定义要对时间变量（这里是times）进行正交多项式转变，我们这里有5个时间点，所以是1次方到4次方：
代码语言：javascript复制contrasts(df.l$times) <- contr.poly(5)
contrasts(df.l$times)
## .L .Q .C ^4
## t0 -6.324555e-01 0.5345225 -3.162278e-01 0.1195229
## t1 -3.162278e-01 -0.2672612 6.324555e-01 -0.4780914
## t2 -3.510833e-17 -0.5345225 1.755417e-16 0.7171372
## t3 3.162278e-01 -0.2672612 -6.324555e-01 -0.4780914
## t4 6.324555e-01 0.5345225 3.162278e-01 0.1195229
然后继续进行方差分析，此时是单纯探索时间对因变量的影响，所以注意formula的形式：
代码语言：javascript复制# A组
f1 <- aov(hp ~ times, data = df.l[df.l$group=="A",])
# 分别看不同次方的结果
summary(f1,
split=list(times=list(liner=1,quadratic=2,cubic=3,biquadrate=4)))
## Df Sum Sq Mean Sq F value Pr(>F)
## times 4 475.4 118.9 5.580 0.003486 **
## times: liner 1 84.5 84.5 3.967 0.060229 .
## times: quadratic 1 26.4 26.4 1.240 0.278655
## times: cubic 1 364.5 364.5 17.113 0.000511 ***
## times: biquadrate 1 0.0 0.0 0.001 0.972627
## Residuals 20 426.0 21.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# B组
f2 <- aov(hp ~ times, data = df.l[df.l$group=="B",])
summary(f2, split=list(times=list(liner=1,quadratic=2,cubic=3,biquadrate=4)))
## Df Sum Sq Mean Sq F value Pr(>F)
## times 4 1017.0 254.3 9.757 0.000152 ***
## times: liner 1 662.5 662.5 25.421 6.24e-05 ***
## times: quadratic 1 296.2 296.2 11.367 0.003034 **
## times: cubic 1 3.9 3.9 0.150 0.702229
## times: biquadrate 1 54.4 54.4 2.088 0.163954
## Residuals 20 521.2 26.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# C组
f3 <- aov(hp ~ times, data = df.l[df.l$group=="C",])
summary(f3, split=list(times=list(liner=1,quadratic=2,cubic=3,biquadrate=4)))
## times
## Df Sum Sq Mean Sq F value Pr(>F)
## times 4 1681.6 420.4 40.915 3.28e-08 ***
## times: liner 1 403.3 403.3 39.249 1.13e-05 ***
## times: quadratic 1 41.7 41.7 4.054 0.0612 .
## times: cubic 1 605.5 605.5 58.931 9.43e-07 ***
## times: biquadrate 1 631.1 631.1 61.425 7.23e-07 ***
## Residuals 16 164.4 10.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以看到结果和课本差异很大，关于这方面的资料较少，如果有朋友知道，欢迎指教！
时间点比较
课本说因为事后检验重复次数太多难以承受，但是我们用计算机很快，所以用事后检验也没什么问题。
事后检验可以参考组间比较，根据组别进行分组，分组比较不同时间点的差别。
事前检验课本采用配对t检验，全都和t0的数据进行比较。
事前检验使用rstatix包解决:
library(rstatix)
df.l |>
group_by(group) |>
t_test(hp ~ times, ref.group = "t0",paired = T)
## # A tibble: 12 × 11
## group .y. group1 group2 n1 n2 statistic df p p.adj p.adj…¹
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 A hp t0 t1 5 5 8.35 4 1 e-3 4 e-3 **
## 2 A hp t0 t2 5 5 1.77 4 1.52e-1 3.04e-1 ns
## 3 A hp t0 t3 5 5 -3.64 4 2.2 e-2 6.6 e-2 ns
## 4 A hp t0 t4 5 5 0.147 4 8.9 e-1 8.9 e-1 ns
## 5 B hp t0 t1 5 5 1.72 4 1.6 e-1 1.6 e-1 ns
## 6 B hp t0 t2 5 5 4.35 4 1.2 e-2 2.4 e-2 *
## 7 B hp t0 t3 5 5 -8.37 4 1 e-3 3 e-3 **
## 8 B hp t0 t4 5 5 -16.7 4 7.47e-5 2.99e-4 ***
## 9 C hp t0 t1 5 5 1.44 4 2.23e-1 2.92e-1 ns
## 10 C hp t0 t2 5 5 4.75 4 9 e-3 2.8 e-2 *
## 11 C hp t0 t3 5 5 -5.12 4 7 e-3 2.8 e-2 *
## 12 C hp t0 t4 5 5 -1.80 4 1.46e-1 2.92e-1 ns
## # … with abbreviated variable name ¹p.adj.signif
直接给出3组的结果，和课本一模一样~
