使用R做方差分析实现多重比较可视化结果

说明：本文章中为作者R学习笔记，资料及操作流程均来源网络，侵权删！

本文源代码“使用R做方差分析源代码.R”

1. 方差分析假定：正态性（否则建立广义线性模型），独立性（否则建立混合线性模型，定义G矩阵和R矩阵），齐次性（否则混合线性模型，定义G矩阵和R矩阵）

2. 单因素方差分析（为什么高级心统老师讲“边际均值比较”更常用？）

2.1 安装相关R包，并找出数据（来源“agridat，将数据命名dat）。这里使用devtools下载github上的文件，devtools后面的格式是install_github("Package Name","Author Name")，::的作用可以理解成当有多个包下有同 一名字函数时，可以用::指定包

# install.packages("agridat")
# install.packages("devtools")
# devtools::install_github("kwstat/agridat")#这一步是利用devtools下载github上的数据

library(agridat)

#单因素方差分析-数据
data(lasrosas.corn)
dat <- lasrosas.corn
str(dat)

2.2 方差分析代码格式：“m1 = aov(yield ~ nf, data=dat)”，m1为模型保存的名称，aov为R中的方差分析代码，yield为因变量，~波浪号前后分开Y与X变量，nf为X变量（种类），data=表示数据库用哪一个，dat为数据库名称

m1 = aov(yield ~ nf, data=dat)
summary(m1)

2.3 这里在查看数据库时，前面用到str()函数，它是用来紧凑的显示对象内部结构，也就是对象里面有什么，除此之外还可以用head()函数表示，但这种方式只可以查看数据前6行数据。summary()函数可以获取描述性统计量，也就是显示出方程或者其他分析的结果

3. 单因素区组

代码格式：“m2 = aov(yield ~ block +trt, data=dat”（+分割预测变量）

#单因素随机区组-建模
m2 = aov(yield ~ block +trt, data=dat)
summary(m2)

4. 两因素方差分析（无交互）-用新的数据库

#两因素方差分析（无交互）-数据库
data(lucas.switchback)
dat <-lucas.switchback
str(dat)
#两因素方差分析（无交互）-建模
m3 = aov(yield ~ block +trt +period, data=dat)
summary(m3)

代码格式：“m3 = aov(yield ~ block +trt +period, data=dat)”（这里是将block也就是区组作为一个自变量了？而且这个区组是谁的区组？）

5. 两因素方差分析（有交互）

#两因素方差分析（有交互）-数据库
data(lucas.switchback)
dat <- lucas.switchback
#两因素方差分析（有交互）-建模
m4 = aov(yield ~ block +trt*period, data=dat)
summary(m4)

代码格式：“m4 = aov(yield ~ block +trt*period, data=dat”，这里trt*period是R中公式的简写，表示trt + period + trt: period，其中trt: period表示交互作用（: 表示交互项，而*表示所有可能的交互项，比如说y ~ A*B*C可展开为y ~ A+B+C+A:B+A:C+B:C+A:B:C；还可以用.它表示包含除因变量外的所有变量，如果一个数据库中有A、B、C三列预测变量，那y.表示y ~ A+B+C；还可以用^表示交互项达到次数，如y~(A+B+C)^2表示y~A+B+C+A:B+A:C+B:C（奇怪这里不考虑A:B:C吗））（另外这里注意预测变量放置位置，一般来讲，越基础的效应越需要放在表达式前面，首先是协变量，然后是主效应，然后是交互效应，放置位置不同，可能会对结果产生较大的影响（这里的block也属于协变量吗））

6. 多因素方差分析

#多因素方差分析-数据库
data(lucas.switchback)
dat <- lucas.switchback
#多因素方差分析-建模
m5 = aov(yield ~ block + trt*period +cow, data=dat)
summary(m5)

代码格式：“m5 = aov(yield ~ block + trt*period + cow, data=dat)”（这里面为什么只考虑了trt和period之间的交互呢？）

7. 数据正态性检验（这一步应该在前面，因为只有经历验证后是正态、独立、方差齐性的个体才能够进行方差分析），在这里可以使用球形检验（Shapiro-Wilk）检验，也可以用qqplot查看残差的图（这个不是在说明残差是否符合正态分布吗？），判断数据的正态性，也可以做直方图查看

#正态性检验/齐性检验-数据库
data(npk)
dat <- npk
str(dat)

7.1 做直方图：（一般分析时，仅对Y变量进行正态性检验，如果是单因素或者多因素，也可以根据因素分组进行正态性检验；数据量大时，对于稍微偏态的数据，即使不太符合正态分布，也不影响结论）

#正态性检验-画直方图
hist(dat$yield)

7.2 做qq图：（使用car软件包中的qqplot函数）“library(car)    qqPlot(dat$yield)”（qqPlot的大小写一定要写对，否则会出错），这里可以看到数据是否落在置信区间里面

#正态性检验-画qq图
library(car)
qqPlot(dat$yield)
shapiro.test(dat$yield)

7.3 使用Shapiro-wilk检验数据正态分布：“shapiro.test(dat$yield)”，可以根据结果中的P值；来看，如果大于0.05，则报名数据符合正态分布

#使用Shapiro-wilk检验数据正态分布
shapiro.test(dat$yield)

8. 齐性检验，可以使用Bartlet检验和Levene检验、

8.1 Bartlett检验（对数据的正态性非常敏感）：“Bartlett.test(yield ~ N, data=dat)”（这里不清楚~ N是什么意思？），可以根据P值判断，如果P大于0.05，则表明方差齐性

#齐性检验-Bartlett检验
bartlett.test(yield ~ N, data=dat)

8.2 Levene检验（非参数检验方法，使用范围更广）：“leveneTest(yield ~ N, data=dat)”

#齐性检验-LeveneTest检验
library(car)
leveneTest(yield ~ N, data = dat)

9. 多重比较

#安装多重比较的包-LSD
install.packages("agricolae")
#多重比较-数据库
data(npk)
dat <- npk
str(dat)

9.1 首先进行方差分析，以对N进行单因素方差分析，block为区组，代码格式：“m9 = aov(yield ~ block + N, data=dat)”，结果显示显著

#对N进行单因素方差分析，block为区组
m9 = aov(yield ~ block + N, data = dat)
summary(m9)

9.2 LSD多重比较，代码格式：“re1 = LSD.test(m9,”N”)    re1$groups”（这里注意，需要下载agricolae包，并用library加载包），结果可以显示N中两个水平对应的yield值（这里注意，本来输出re1会出来很多结果的，但因为只看groups的结果，所以用了代码re1$groups），应该是可以通过这个进行比较，不过不清楚后面的a和b对应于什么关系，以及没有P值呈现，如果是三个变量是不是就不好判断了？

#LSD
library(agricolae)
re1 = LSD.test(m9, "N")

Scheffe多重比较，代码格式：“(scheffe.test(m9,”N”))”，出来的结果比较多，不过最下面的结果和LSD的是一致的

#scheffe多重比较
(scheffe.test(m9, "N"))

Tukey多重比较（保守，控制所有比较的误差），代码格式：“(HSD.test(m9, “N”)”，结果与scheffe的结果比较像

#Tukey多重比较
(HSD.test(m9,"N"))

SNK多重比较（不如Tukey保守，只控制要比较的组），代码格式：“(SNK.test(m9, “N”))”

#SNK多重比较
(SNK.test(m9,"N"))

Duncan多重比较（只比LSD稍微保守一点）：“(Duncan.test(m9,”N”))”（这里只是自己实验了一下::的使用）

#duncan多重比较
aaa = agricolae::duncan.test(m9,"N")
aaa

10. 多重比较可视化：首先，需要对数据进行变化，由宽数据（数据集对所有变量进行了明确的细分，各变量的值不存在重复循环的情况也无法归类，总体表现为变量多而观察值少）变化为长数据（数据集中的变量没有做明确的细分，即变量中至少有一个变量中的元素存在值严重重复循环的情况，这种可以归为几类），表格整体的形状为长方形，即变量少儿观察值多。变化原因：宽数据无法用ggplot作图。（因为我使用的数据不用转化，所以就没有搞）

10.1 首先找到数据（dat）并建模，得到显著结果后，进行LSD多重比较（这里需要用到agricolae）。

#单因素方差分析-数据
data(lasrosas.corn)
dat <- lasrosas.corn
#单因素方差分析-建模
m1 = aov(yield ~ nf, data=dat)
summary(m1)
#多重比较可视化，dc=多重比较
library(agricolae)
dc = LSD.test(m1,"nf", alpha = 0.05)
dc$groups

10.2 因为作图需要标准误，所以先计算标准误，使用aggregate。

#计算标准误，bw = 标准误
bw = aggregate(yield ~ nf, dat, sd)#dat是数据，sd是要计算的标准误
bw#查看标准误计算结果

10.3 计算完成后得到bw文件，重新命名bw文件

names(bw) = c("nf", "sd")#重新对bw的两列数据进行命名,这里的第二列数据的标题，会对应加上标准误那一行的sd，一定不要搞错
bw#查看改变命名后的数据文件

10.4 下载dplyr和tidyr包，合并数据，合并dc$groups和bw的数据

install.packages("dplyr")
install.packages("tidyr")
library(dplyr)
library(tidyr)

#根据nf合并数据，dc$groups是dc文件中的groups结果，%>%管道
dc2 = dc$groups %>% mutate(nf = rownames(dc$groups)) %>% inner_join(., bw, by = "nf")
#这里很重要，不太理解是怎样将dc$groups中的N1等列数据加上nf名称的

10.5 作图

10.5.1直方图

##作直方图
library("ggplot2")
p1 = dc2 %>% ggplot(aes(x=nf,y=yield)) + geom_col(aes(fill = nf),width=.4)
p1

加标准误（这里的sd要注意，和前面重新命名bw那里的sd是一样）

#加上标准误
p2 = p1 +  geom_errorbar(aes(ymax = yield + sd, ymin = yield - sd),width = .1,size = .5)
p2

加显著性

#加上显著性
p3 = p2 +geom_text(aes(label = groups, y=yield + sd +0.5))
p3

背景设置为白色

#背景设置为白色
p4 = p3 + theme(panel.grid = element_blank(), panel.background = element_rect(color = "black",fill = "transparent"))
p4

增加标签

#增加标签
p5 = p4 + labs(x = "NF水平", y = "yield结果",title = "多重比较可视化")
p5

对横轴nf重新排序，并重新上面的几步

#将顺序按照：N1，N2，N3，N4，N5，N0的顺序
dc2$nf = factor(dc2$nf,levels = c("N1","N2","N3","N4","N5","N0"))
#重新进行上面的几步
p1 = dc2 %>% ggplot(aes(x=nf,y=yield)) + geom_col(aes(fill = nf),width=.4)
p2 = p1 +  geom_errorbar(aes(ymax = yield + sd, ymin = yield - sd),width = .1,size = .5)
p3 = p2 +geom_text(aes(label = groups, y=yield + sd +0.5))
p4 = p3 + theme(panel.grid = element_blank(), panel.background = element_rect(color = "black",fill = "transparent"))
p5 = p4 + labs(x = "NF水平", y = "yield结果",title = "多重比较可视化")
p5

加趋势线

#加趋势线
p6 = p5 +geom_line(aes(group=""),color = "red")
p6

补充知识点

1.LSD本质是做了多次t检验，同时没有进行矫正，所以一般不用

2.当有协变量时，事后检验不能做

3.事后多重比较的t值的df是方差分析表中error的df（有的报告t有的报告F，F其实是t的平方，它的自由度还是error中的df

4.先验比较（操作中点击“对比(N)”）不需要事后检验