今天继续学习使用R语言进行医学统计学分析，今天要学习的内容是协方差分析，还是使用课本中的例子。
我使用的课本是孙振球主编的《医学统计学》第4版，封面如下：
课本封面
完全随机设计资料的协方差分析
使用课本例13-1的例子。
首先是读取数据，本次数据手动录入：
代码语言：javascript复制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)
)
看一下数据结构：
代码语言：javascript复制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 ...
可以看到一共6列，和课本上面的一模一样，分别是x1,y1,x2,y2,x3,y3。
接下来为了进行方差分析，需要变为长数据，把所有的x放在1列，所有的y放在1列，还有一列是组别：
如果大家还对长宽数据转换不了解的，赶紧翻看之前的推文：
长数据变为宽数据的7种情况！
宽数据变为长数据的5种情况！
这是一个非常重要且使用频率极高的技能！
代码语言：javascript复制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 <fct> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1…
## $ x <dbl> 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 <dbl> 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…
所有的x放在1列，所有的y放在1列，还有一列是组别！
然后就是进行单因素协方差分析：
代码语言：javascript复制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
得到的结果和课本是一模一样的，组内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
一行代码即可：
代码语言：javascript复制ancovaplot(y ~ x group, data = df13_11)
plot of chunk unnamed-chunk-6
但其实我们也可以用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'
plot of chunk unnamed-chunk-7
好看是好看，但是很明显不如HH简洁啊！
使用rstatix进行优雅的协方差分析
代码语言：javascript复制library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:stats':
##
## filter
进行协方差分析：
代码语言：javascript复制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
结果也是一样的！
