|
学术论文和杂志中最常见的是用“字母”和“星号”标记差异显著性，在之前的文章《显著性分析后如何标记“abc”？》已经介绍过如何用SPSS比较均值的结果进行“纯手动”添加字母显著性标记。嗯，方法似乎非常“绕”，且一走神就容易错。
那么有没有R脚本可以实现相关的方差分析、均值比较、生成变量均值对应的“字母标记”？当然，如果能直接将字母标记添加到图表(比如柱状图)中，岂不美滋滋？
闲言少叙，接下来主要为大家介绍如何用R进行方差齐性检验（Bartlett test 和Levene test）、方差分析、均值的多重比较方法（TukeyHSD和LSD法），最后用ggplot2包进行数据可视化。示例数据和脚本可通过点击 阅读原文 下载
# 读取示例数据 data<>'1233.txt',sep='\t',header= TRUE) data 方差齐性检验 nom<-bartlett.test(dia~lab,data ==""> # 如果是2（多）因数，使用interaction()函数，leveneTest(y~interaction(var1,var2),data= data) nom nom$p.value # 输出p值，当p>0.05时，方差是齐性的。 # 使用另一种方法进行齐性检验，SPSS软件使用Levene test进行齐性检验，这里需要安装car包。
install.packages('car') library('car') nom1<-levenetest(dia~lab,data ==""> # 如果是2因素，使用*，leveneTest(y~ivar1*var2,data = data) nom1 nom1$`Pr(>F)` # 输出p值，当p>0.05时，方差是齐性的。
方差分析 # 单因素方差分析，整体来看差异显著 oneway<-aov(dia~lab,data ==""> anova(oneway)
多重比较 #TukeyHSD法（Tukey’s Honestly Significant Difference） bijiao<-tukeyhsd(oneway,ordered =false,conf.level=""> bijiao # 得到两两比较的差异分析结果
# LSD法（Fisher’s Least Significant Difference） install.packages('agricolae') library('agricolae') out <->'lab',p.adj='bonferroni') out # LSD法检验处微小的差异，比较方便的是直接得出显著行标记，不需人工标记 # 整理用于作图的数据框
mar<> rownamemar<> newmar<> sort<> # 将groups的数据框按列名排序，目的是保持与均值标准差的数据一一对应
rowname<> mean<> sd<> marker<> plotdata<> plotdata # 生成新的数据框
作图+显著性标记 # ggplot2 绘制带显著性标记的柱状图，更详细的内容参考《ggplot2绘图学习笔记分享》的内容
library('ggplot2') p1<-ggplot(plotdata,aes(x=factor(rowname),y=mean))+geom_bar(position=position_dodge(0.6),width =="" 0.5,stat="">'identity') p1 p2<> p3<-p2+geom_text(aes(x=factor(rowname),y=mean+sd+2.0,label=marker),size=3,position=> p3 p4<>'Lesiondiameter (mm)') p4 p5<-p4+coord_cartesian(ylim= c(0,60),expand=""> p5
#更改y轴显示范围，这里的expand默认为TRUE mytheme<-theme_bw()+theme(axis.title =element_text(size=""> axis.text =element_text(size=12), panel.grid.major =element_line(color = 'white'), panel.grid.minor =element_line(colour = 'white'), axis.text.x =element_text(size = 12,angle=0,vjust=0,hjust=0,color = 'black'), axis.text.y =element_text(size = 12,color = 'black'), legend.text =element_text(size = 12),legend.title = element_blank() plot.margin =unit(c(0.5,0.5,0.5,0.5), 'cm')) p5+mytheme
# 最终的绘制结果如下，如果不喜欢柱子下的“空隙”，可将p5中coord_cartesian(ylim = c(0,60),expand =FALSE)的expand改为TRUE。 #将图表保存为PDF文件 ggsave('wanyixia1.pdf', width = 10, height= 10, units = 'cm')
最后想一下，如果不用LSD法现成的标记，用TukeyHSD法两两比较的结果又该进行标记呢？
今天的内容就到这里啦~ |
|
|
