R语言和医学统计学（13）：协方差分析

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

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

今天继续学习使用R语言进行医学统计学分析，今天要学习的内容是协方差分析，还是使用课本中的例子。

我使用的课本是孙振球主编的《医学统计学》第4版，封面如下：

完全随机设计资料的协方差分析

使用课本例13-1的例子。

首先是读取数据，本次数据手动录入：

df13_1 <- data.frame(x1=c(10.8,11.6,10.6,9.0,11.2,9.9,10.6,10.4,9.6,10.5,10.6,9.9,9.5,9.7,10.7,9.2,10.5,11.0,10.1,10.7,8.5,10.0,10.4,9.7,9.4,9.2,10.5,11.2,9.6,8.0),
                     y1=c(9.4,9.7,8.7,7.2,10.0,8.5,8.3,8.1,8.5,9.1,9.2,8.4,7.6,7.9,8.8,7.4,8.6,9.2,8.0,8.5,7.3,8.3,8.6,8.7,7.6,8.0,8.8,9.5,8.2,7.2),
                     x2=c(10.4,9.7,9.9,9.8,11.1,8.2,8.8,10.0,9.0,9.4,8.9,10.3,9.3,9.2,10.9,9.2,9.2,10.4,11.2,11.1,11.0,8.6,9.3,10.3,10.3,9.8,10.5,10.7,10.4,9.4),
                     y2=c(9.2,9.1,8.9,8.6,9.9,7.1,7.8,7.9,8.0,9.0,7.9,8.9,8.9,8.1,10.2,8.5,9.0,8.9,9.8,10.1,8.5,8.1,8.6,8.9,9.6,8.1,9.9,9.3,8.7,8.7),
                     x3=c(9.8,11.2,10.7,9.6,10.1,9.8,10.1,10.3,11.0,10.5,9.2,10.1,10.4,10.0,8.4,10.1,9.3,10.5,11.1,10.5,9.7,9.2,9.3,10.4,10.0,10.3,9.9,9.4,8.3,9.2),
                     y3=c(7.6,7.9,9.0,7.8,8.5,7.5,8.3,8.2,8.4,8.1,7.0,7.7,8.0,6.6,6.1,8.1,7.8,8.4,8.2,8.0,7.6,6.9,6.7,8.1,7.4,8.2,7.6,7.8,6.6,7.2)
                     )
1
2
3
4
5
6
7

看一下数据结构：

str(df13_1)
## 'data.frame':	30 obs. of  6 variables:
##  $ x1: num  10.8 11.6 10.6 9 11.2 9.9 10.6 10.4 9.6 10.5 ...
##  $ y1: num  9.4 9.7 8.7 7.2 10 8.5 8.3 8.1 8.5 9.1 ...
##  $ x2: num  10.4 9.7 9.9 9.8 11.1 8.2 8.8 10 9 9.4 ...
##  $ y2: num  9.2 9.1 8.9 8.6 9.9 7.1 7.8 7.9 8 9 ...
##  $ x3: num  9.8 11.2 10.7 9.6 10.1 9.8 10.1 10.3 11 10.5 ...
##  $ y3: num  7.6 7.9 9 7.8 8.5 7.5 8.3 8.2 8.4 8.1 ...
1
2
3
4
5
6
7
8

可以看到一共6列，和课本上面的一模一样，分别是x1,y1,x2,y2,x3,y3。

接下来为了进行方差分析，需要变为长数据，把所有的x放在1列，所有的y放在1列，还有一列是组别：

如果大家还对长款数据转换不了解的，赶紧翻看之前的推文：

dsfasdfasdfaf

这是一个非常重要且使用频率极高的技能！

suppressPackageStartupMessages(library(tidyverse))

df13_11 <- df13_1 %>% 
  pivot_longer(cols = everything(), # 变长
               names_to = c(".value","group"),
               names_pattern = "(.)(.)"
               ) %>% 
  mutate(group = as.factor(group)) # 组别变为因子型

glimpse(df13_11) # 查看数据结构，神奇！
## Rows: 90
## Columns: 3
## $ group  1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1…
## $ x      10.8, 10.4, 9.8, 11.6, 9.7, 11.2, 10.6, 9.9, 10.7, 9.0, 9.8, 9.6…
## $ y      9.4, 9.2, 7.6, 9.7, 9.1, 7.9, 8.7, 8.9, 9.0, 7.2, 8.6, 7.8, 10.0…
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

所有的x放在1列，所有的y放在1列，还有一列是组别！

然后就是进行单因素协方差分析：

fit <- aov(y ~ x + group, data = df13_11) # 注意公式的写法，一定是把协变量放在主变量前面！
summary(fit)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## x            1  29.06  29.057  171.20 <2e-16 ***
## group        2  19.85   9.925   58.48 <2e-16 ***
## Residuals   86  14.60   0.170                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1
2
3
4
5
6
7
8

得到的结果和课本是一模一样的，组内ss=14.60, ms=0.170, v=86, 修正均数ss=19.85, ms=9.925, v=2, F=58.48，拒绝H0，接受H1，可以认为在扣除初始（基线）糖化血红蛋白含量的影响后，3组患者的总体降糖均数有差别。

结果的可视化可以使用HH包：

library(HH)
## Loading required package: lattice
## Loading required package: grid
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'HH'
## The following object is masked from 'package:purrr':
## 
##     transpose
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
35

一行代码即可：

ancovaplot(y ~ x + group, data = df13_11)
1

但其实我们也可以用ggplot2来画，可能更好看一点：

theme_set(theme_minimal())

p1 <- ggplot(df13_11, aes(x=x,y=y))+
  geom_point(aes(color=group,shape=group))+
  geom_smooth(method = "lm",se=F,aes(color=group))+
  labs(y=NULL)

p2 <- ggplot(df13_11, aes(x=x,y=y))+
  geom_point(aes(color=group,shape=group))+
  geom_smooth(method = "lm",se=F,aes(color=group))+
  facet_wrap(~group)

library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
p2 + p1 + plot_layout(guides = 'collect',widths = c(3, 1))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

好看是好看，但是很明显不如HH简洁啊！

使用rstatix进行优雅的协方差分析

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
1
2
3
4
5
6
7
8
9

进行协方差分析：

res <- anova_test(y ~ x + group, data = df13_11, type = 1)
get_anova_table(res)
## ANOVA Table (type I tests)
## 
##   Effect DFn DFd       F        p p<.05   ges
## 1      x   1  86 171.199 3.64e-22     * 0.666
## 2  group   2  86  58.480 9.22e-17     * 0.576
1
2
3
4
5
6
7

结果也是一样的！

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

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