R语言标记显著性标记

学术论文和杂志中最常见的是用“字母”和“星号”标记差异显著性,在之前的文章《显著性分析后如何标记“abc”?》已经介绍过如何用SPSS比较均值的结果进行“纯手动”添加字母显著性标记。嗯,方法似乎非常“绕”,且一走神就容易错。

那么有没有R脚本可以实现相关的方差分析、均值比较、生成变量均值对应的“字母标记”?当然,如果能直接将字母标记添加到图表(比如柱状图)中,岂不美滋滋?

闲言少叙,接下来主要为大家介绍如何用R进行方差齐性检验(Bartlett test 和Levene test)、方差分析、均值的多重比较方法(TukeyHSD和LSD法),最后用ggplot2包进行数据可视化。示例数据和脚本可通过点击 阅读原文 下载

# 读取示例数据

  1. data<-read.table("1233.txt",sep="t",header= TRUE)
  2. data

R语言标记显著性标记-图片1

方差齐性检验

  1. nom<-bartlett.test(dia~lab,data = data)
  2. # 如果是2(多)因数,使用interaction()函数,leveneTest(y~interaction(var1,var2),data= data)
  3. nom
  4. nom$p.value # 输出p值,当p>0.05时,方差是齐性的。

R语言标记显著性标记-图片2

  1. # 使用另一种方法进行齐性检验,SPSS软件使用Levene test进行齐性检验,这里需要安装car包。
  2. install.packages("car")
  3. library("car")
  4. nom1<-leveneTest(dia~lab,data = data)
  5. # 如果是2因素,使用*,leveneTest(y~ivar1*var2,data = data)
  6. nom1
  7. nom1$`Pr(>F)` # 输出p值,当p>0.05时,方差是齐性的。

R语言标记显著性标记-图片3

方差分析

  1. # 单因素方差分析,整体来看差异显著
  2. oneway<-aov(dia~lab,data = data)
  3. anova(oneway)

R语言标记显著性标记-图片4

多重比较

  1. #TukeyHSD法(Tukey’s Honestly Significant Difference)
  2. bijiao<-TukeyHSD(oneway,ordered =FALSE,conf.level = 0.95)
  3. bijiao
  4. # 得到两两比较的差异分析结果

R语言标记显著性标记-图片5

  1. # LSD法(Fisher’s Least Significant Difference)
  2. install.packages("agricolae")
  3. library("agricolae")
  4. out <- LSD.test(oneway,"lab",p.adj="bonferroni")
  5. out
  6. # LSD法检验处微小的差异,比较方便的是直接得出显著行标记,不需人工标记

R语言标记显著性标记-图片6

  1. # 整理用于作图的数据框
  2. mar<-out$groups
  3. rownamemar<-row.names(mar)
  4. newmar<-data.frame(rownamemar,mar$dia,mar$groups)
  5. sort<-newmar[order(newmar$rownamemar),]
  6. # 将groups的数据框按列名排序,目的是保持与均值标准差的数据一一对应
  7. rowname<-row.names(out$means)
  8. mean<-out$means[,1]
  9. sd<-out$means[,2]
  10. marker<-sort$mar.groups
  11. plotdata<-data.frame(rowname,mean,sd,marker)
  12. plotdata
  13. # 生成新的数据框

R语言标记显著性标记-图片7

作图+显著性标记

  1. # ggplot2 绘制带显著性标记的柱状图,更详细的内容参考《ggplot2绘图学习笔记分享》的内容
  2. library("ggplot2")
  3. p1<-ggplot(plotdata,aes(x=factor(rowname),y=mean))+geom_bar(position=position_dodge(0.6),width = 0.5,stat = "identity")
  4. p1
  5. p2<-p1+geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),position=position_dodge(0.6),width=0.2)
  6. p3<-p2+geom_text(aes(x=factor(rowname),y=mean+sd+2.0,label=marker),size=3,position= position_dodge(0.6))
  7. p3
  8. p4<-p3+xlab("")+ylab("Lesiondiameter (mm)")
  9. p4
  10. p5<-p4+coord_cartesian(ylim= c(0,60),expand = FALSE)
  11. p5
  12. #更改y轴显示范围,这里的expand默认为TRUE
  13. mytheme<-theme_bw()+theme(axis.title =element_text(size = 12),
  14. axis.text =element_text(size=12),
  15. panel.grid.major =element_line(color ="white"),
  16. panel.grid.minor =element_line(colour = "white"),
  17. axis.text.x =element_text(size = 12,angle=0,vjust=0,hjust=0,color = "black"),
  18. axis.text.y =element_text(size = 12,color ="black"),
  19. legend.text =element_text(size = 12),legend.title = element_blank()
  20. plot.margin =unit(c(0.5,0.5,0.5,0.5), "cm"))
  21. p5+mytheme
  22. # 最终的绘制结果如下,如果不喜欢柱子下的“空隙”,可将p5中coord_cartesian(ylim = c(0,60),expand =FALSE)的expand改为TRUE。

 

R语言标记显著性标记-图片8

  1. #将图表保存为PDF文件
  2. ggsave("wanyixia1.pdf", width = 10, height= 10, units = "cm")

最后想一下,如果不用LSD法现成的标记,用TukeyHSD法两两比较的结果又该进行标记呢?

发表评论

匿名网友

拖动滑块以完成验证
加载失败