R语言和医学统计学（6）：重复测量方差分析

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

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

前言

这是R语言和医学统计学的第6篇内容。

主要是用R语言复现课本中的例子。我使用的课本是孙振球主编的《医学统计学》第4版，封面如下：

重复测量数据两因素两水平的方差分析

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

df12_1 <- foreign::read.spss("E:/各科资料/医学统计学/研究生课程/析因设计重复测量/9重复测量18-9研/12-1.sav", to.data.frame = T)

str(df12_1)
## 'data.frame':	20 obs. of  5 variables:
##  $ n    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ x1   : num  130 124 136 128 122 118 116 138 126 124 ...
##  $ x2   : num  114 110 126 116 102 100 98 122 108 106 ...
##  $ group: Factor w/ 2 levels "处理组","对照组": 1 1 1 1 1 1 1 1 1 1 ...
##  $ d    : num  16 14 10 12 20 18 18 16 18 18 ...
##  - attr(*, "variable.labels")= Named chr [1:5] "编号" "治疗前血压" "治疗后血压" "组别" ...
##   ..- attr(*, "names")= chr [1:5] "n" "x1" "x2" "group" ...
##  - attr(*, "codepage")= int 936
1
2
3
4
5
6
7
8
9
10
11
12

数据一共5列（第5列是自己算出来的，其实原始数据只有4列），第1 列是编号，第2列是治疗前血压，第3例是治疗后血压，第4列是分组，第5列是血压前后差值。

进行重复测量数据两因素两水平的方差分析前，先把数据转换一下格式：

library(tidyverse)
## -- Attaching packages ----------------------------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.3     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## -- Conflicts -------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
df12_11 <- 
  df12_1[,1:4] %>% 
  pivot_longer(cols = 2:3,names_to = "time",values_to = "hp") %>% 
  mutate_if(is.character, as.factor)

df12_11$n <- factor(df12_11$n)

str(df12_11)
## tibble [40 x 4] (S3: tbl_df/tbl/data.frame)
##  $ n    : Factor w/ 20 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $ group: Factor w/ 2 levels "处理组","对照组": 1 1 1 1 1 1 1 1 1 1 ...
##  $ time : Factor w/ 2 levels "x1","x2": 1 2 1 2 1 2 1 2 1 2 ...
##  $ hp   : num [1:40] 130 114 124 110 136 126 128 116 122 102 ...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

转换后的数据格式如上，只截取了一部分。

进行重复测量数据两因素两水平的方差分析:

hp是因变量，time是测量时间（治疗前和治疗后各测量一次），group是分组因素（两种治疗方法），n是受试者编号。

f1 <- aov(hp ~ time * group + Error(n/(time)), data = df12_11)

summary(f1)
## 
## Error: n
##           Df Sum Sq Mean Sq F value Pr(>F)
## group      1  202.5   202.5   1.574  0.226
## Residuals 18 2315.4   128.6               
## 
## Error: n:time
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## time        1 1020.1  1020.1   55.01 7.08e-07 ***
## time:group  1  348.1   348.1   18.77 0.000401 ***
## Residuals  18  333.8    18.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

结果输出了两张表，第二个是测量前后比较与交互作用的方差分析表，第一个是处理组与对照组比较的方差分析表，可以看到结果和课本是一样的！

用图形方式展示重复测量的结果：

with(df12_11,
     interaction.plot(time, group, hp, type = "b", col = c("red","blue"), 
                      pch = c(12,16), main = "两因素两水平重复测量方差分析"))
1
2
3

或者用箱线图展示结果：

boxplot(hp ~ group*time, data = df12_11, col = c("gold","green"),
        main = "两因素两水平重复测量方差分析")
1
2

重复测量数据两因素多水平的分析

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

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] "序号" "组别" "" "" ...
##   ..- 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(tidyverse)

df12_31 <- df12_3 %>% 
  pivot_longer(cols = 3:7, names_to = "times", values_to = "hp")

df12_31$No <- factor(df12_31$No)
df12_31$times <- factor(df12_31$times)

str(df12_31)
## tibble [75 x 4] (S3: tbl_df/tbl/data.frame)
##  $ No   : Factor w/ 15 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ group: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
##  $ times: Factor w/ 5 levels "t0","t1","t2",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ hp   : num [1:75] 120 108 112 120 117 118 109 115 126 123 ...
1
2
3
4
5
6
7
8
9
10
11
12
13
14

转换后的格式见上图，只截取了部分。

进行方差分析：

f2 <- aov(hp ~ times * group + Error(No/(times)), data = df12_31)

summary(f2)
## 
## 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 ***
## times:group  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

输出结果是两张表格，第1个是不同诱导方法患者血压比较的方差分析表，第2个是麻醉诱导时相及其与诱导方法交互作用的方差分析表。

结果和课本是一样的哟！具体意义解读请认真学习医学统计学相关知识。

用图形方式展示重复测量的结果：

with(df12_31,
     interaction.plot(times, group, hp, type = "b", 
                      col = c("red","blue","green"), 
                      pch = c(12,16,20), 
                      main = "两因素多水平重复测量方差分析"))
1
2
3
4
5

或者用箱线图展示结果：

boxplot(hp ~ group*times, data = df12_31, col = c("gold","green","black"),
        main = "两因素多水平重复测量方差分析")
1
2

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

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