(一)使用函数对数据进行参数或非参数检验
(1)导入数据,进行正态性检验,方差齐性检验
library(carData)library(car)a<-read.csv("C:/Users/baixin/Desktop/cml2.csv",fileEncoding = "GBK")#表格编码为GBK格式,否则无法写入,rstudio使用UTF-8编码#View(a)#cyl<- subset(a, species == 'species列名称' & treatment %in% c('CK', 'DR2'))#仅选取两组的代码#选取species列名称相同物种的所有处理cyl2<- subset(a, species == 'species列名称')cyl2#多组间正态分布检验cyl2_shapiro <- tapply(cyl2$mfl2, cyl2$treatment, shapiro.test)cyl2_shapiro#结果部分不符合正态分布##Shapiro-Wilk 检验,当且仅当两者 p 值均大于 0.05 时表明数据符合正态分布#检验方差齐性#Levene检验方差齐性,非正态分布与正态分布数据均适用cyl2_treatment_factor<-as.factor(cyl2$treatment)#转为factor类型cyl2_treatment_factor<-factor(cyl2_treatment_factor,levels= c('CK' ,'CK2','DR2',#CK CK10 CK2 CK3 CK8 CK9 DR10 DR2 DR3 DR8 DR9'CK3', 'DR3','CK8','DR8','CK9', 'DR9', 'CK10', 'DR10'))#将因子水平重新排序cyl2_treatment_factorcyl2_leveneTest<-leveneTest(cyl2$mfl2, cyl2_treatment_factor,center=median)#检验方差齐性cyl2_leveneTest
(2)数据部分不符合正态分布,使用非参数检验
#####由于部分数据不符合正态分布,使用非参数检验多重比较install.packages("PMCMRplus")library(PMCMRplus)#检验总体kruskal.test(cyl2$mfl2~cyl2$treatment)#结果表明总体有差异#进行两组间检验cyl2$treatment<-factor(cyl2$treatment,levels= c('CK' ,'CK2','DR2',#CK CK10 CK2 CK3 CK8 CK9 DR10 DR2 DR3 DR8 DR9'CK3', 'DR3','CK8','DR8','CK9', 'DR9' ,'CK10', 'DR10'))#将因子水平重新排序feican_test<-kwAllPairsDunnTest(cyl2$mfl2~cyl2$treatment,cyl2,p.adjust.method = "bonferroni")feican_test#会有警告信息##全部组进行非参数检验导致全无差异,可以将数据分为单个比较对,分别使用参数和非参数检验进行尝试?
(3)数据符合正态分布,则进行t.test参数检验。
###独立样本的 t 检验group_name <- unique(cyl2$treatment)number<-length(group_name)#计算共分为几组t_result_two <- NULLfor (i in 1:(number - 1)) {#从范围中选择i和jfor (j in (i + 1):number) {group_ij <- subset(cyl2, treatment %in% c(group_name[i], group_name[j]))#选择两个分组t_result_otu_ij <- t.test(mfl2~treatment, data = group_ij, paired =FALSE ,alternative = 'two.sided',var.equal = TRUE)#t检验t_result_two <- rbind(t_result_two,c(paste(group_name[i], group_name[j], sep = '/'),t_result_otu_ij$p.value))}}#var.equal = FALSE方差不齐,alternative = 'two.sided'两侧检验#paired = FALSE两样本独立t_result_two <- data.frame(t_result_two,check.names = FALSE, stringsAsFactors = F)#t_result_twonames(t_result_two)<-c("两样本","p-value")#添加列名t_result_twofor (i in 1:nrow(t_result_two)) {if (t_result_two[i, 'p-value'] <= 0.001) t_result_two[i, 'SIG'] <- '***'else if (t_result_two[i, 'p-value'] <= 0.01) t_result_two[i, 'SIG'] <- '**'else if (t_result_two[i, 'p-value'] <= 0.05) t_result_two[i, 'SIG'] <- '*'}t_result_twowrite.csv(t_result_two ,"C:/Users/baixin/Desktop/t_result_two.csv")#将结果写出到csv文件
(二)直接在ggplot中进行绘制箱线图,并进行显著性检验
####ggplot绘制带有显著性标记的箱线图,直接能画出非参数检验的结果,不用单独检验,比较要用到ggpubr包。library(ggpubr)library(ggplot2)cyl2$treatment<-factor(cyl2$treatment,levels= c('CK' ,'CK2','DR2',#CK CK10 CK2 CK3 CK8 CK9 DR10 DR2 DR3 DR8 DR9'CK3', 'DR3','CK8','DR8','CK9', 'DR9' ,'CK10', 'DR10'))#将因子水平重新排序my_comparisons <- list( c("DR2", "CK2"),#构建对比矩阵c("DR3", "CK3"),c("DR8", "CK8"),c("DR9", "CK9"),c("DR10", "CK10"))box_plot_per<-ggplot(data=cyl2, aes(x=treatment,y=mfl,fill=treatment,color=treatment) )box_plot_perbox_plot<-box_plot_per+geom_boxplot()+geom_jitter(width = 0.3, size=0.01)+stat_compare_means(method = "kruskal.test", comparisons = my_comparisons,label = "p.signif")#method表示所用检验模型,kruskal.test为非参数检验,t.test为参数检验box_plotbox_plot_1<-box_plot+theme_bw() +#去除背景色theme(panel.grid=element_blank())#去除网格线###ggplot绘制结束
(三)也可以画出分组箱线图,再将对比结果进行添加
#画随时间变化的趋势图p_compare <-ggplot(cyl2, aes(x = time, y = mfl2,fill=group_treatment))+#aes框内fill以fill分组控制填充颜色,color控制线颜色,aes框外无法进行分组geom_boxplot(aes())+theme_bw()+ #去除背景色theme(panel.grid=element_blank())p_compare
结果如图