R语言重复测量方差分析的多重比较

---

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

---

前面介绍了多个样本均数的多重比较，多样本非参数检验后的多重比较。

今天学习下重复测量数据的多重比较，本篇内容和课本结果差异较大，如有错误欢迎指出。

使用的课本是孙振球，徐勇勇《医学统计学》第4版。

重复测量方差分析

使用课本例12-1的数据，直接读取：

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] "\xd0\xf2\xba\xc5" "\xd7\xe9\xb1\xf0" "" "" ...
##   ..- attr(*, "names")= chr [1:7] "No" "group" "t0" "t1" ...
1
2
3
4
5
6
7
8
9
10
11
12
13

数据一共7列，第1列是患者编号，第2列是诱导方法（3种），第3-7列是5个时间点的血压。

首先转换数据格式：

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

进行重复测量方差分析，默认方法不能输出球形检验的结果，所以我更推荐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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 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         *
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

画图展示：

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()
1
2
3
4
5
6
7
8

接下来是重复测量数据的多重比较，课本中分成了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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

P值和课本不太一样，但是结论是一样的，A组和B组之间，A组和C组之间有差别，B组和C组之间没有差别。

时间趋势比较

重复测量方差分析可以采取正交多项式来探索时间变化趋势，具体的内涵解读可以参考冯国双老师的这篇文章：https://mp.weixin.qq.com/s/ndinwbDJsHjAelvNfwqgwA

在R里面进行正交多项式的探索略显复杂，首先定义要对时间变量（这里是times）进行正交多项式转变，我们这里有5个时间点，所以是1次方到4次方：

contrasts(df.l$times) <- contr.poly(5,scores = c(0,1,2,3,4))
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
1
2
3
4
5
6
7
8

然后继续进行方差分析，此时是单纯探索时间对因变量的影响，所以注意formula的形式：

# 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 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
1
2
3
4
5
6
7
8
9
10
11
12

# C组
f3 <- aov(hp ~ times+Error(No/times), data = df.l[df.l$group=="C",])
summary(f3, split=list(times=list(liner=1,quadratic=2,cubic=3,biquadrate=4)))
## 
## Error: No
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  4     98    24.5               
## 
## Error: No: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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

可以看到结果和课本差异很大，关于这方面的资料较少，如果有朋友知道，欢迎指教！

时间点比较

课本说因为事后检验重复次数太多难以承受，但是我们用计算机很快，所以用事后检验也没什么问题。

事后检验可以参考组间比较，根据组别进行分组，分组比较不同时间点的差别。

事前检验课本采用配对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…¹
##  *                       
##  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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

直接给出3组的结果，和课本一模一样~

---

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

---