R语言基础数据分析—单因素方差分析

有了试验数据，我们就需要进行数据的处理与分析，而在试验设计中，通常分为单因素试验或者双因素试验。试验中要考察的指标称为试验指标，影响试验指标的条件称为因素，因素所处的状态称为水平，若试验中只有一个因素改变则称为单因素试验，若有两个因素改变则称为双因素试验，若有多个因素改变则称为多因素试验。方差分析就是对试验数据进行分析，检验方差相等的多个正态总体均值是否相等，进而判断各因素对试验指标的影响是否显著，根据影响试验指标条件的个数可以区分为单因素方差分析、双因素方差分析和多因素方差分析。

在进行数据分析之前，我们要对试验数据进行检验，即独立性检验、正态性检验、方差齐性检验，若三者检验符合条件，则我们可以进行参数分析，如果不符合条件则只能进行非参数分析。

（参数分析：假定数据服从某分布（一般为正态分布），通过样本参数的估计量（x±s）对总体参数（μ）进行检验，比如t检验、方差分析、Pearson的相关分析等。非参数分析：不需要假定总体分布形式，直接对数据的分布进行检验。由于不涉及总体分布的参数，故名「非参数」检验，包括卡方检验、二项分布检验、K-S检验、秩和检验、符号检验、Spearman和Kendall的相关性分析等）

接下来，我演示一个R语言进行的参数分析的单因素方差分析。

原始数据：探究添加不同浓度的氮肥（N1,N2,N3,N4,N5）对植物根际土壤硝态氮的浓度影响

（同学们使用edit（）函数，把我的数据仿照输入，这样你们的数据不需要整理，直接进入数据重构的那步。）

data<-data.frame(N1=numeric(0),N2=numeric(0),N3=numeric(0),N4=numeric(0),N5=numeric(0))
data<-edit(data)

第一步：进行数据整理与重构：

library(readxl)
DN <- read_excel("C:/Users/lenovo/Desktop/DN.xlsx")
View(DN)
DN
data<-edit(DN)#打开数据编辑表，查看和编辑数据

发现数据错误，进行数据的调整（见上一文章）

library(tidyverse)
library(gcookbook)
data<-DN%>%select(N1,N2,N3,N4,N5)
data

数据重构：

library(reshape2)
data <- melt(data,id=c())#进行数据的合并重构
names(data) <- c('trt', 'val') 
data

第二步：对数据进行参数检验（独立性、正态性和方差齐性检验）

data1<-table(data$trt,data$val)
chisq.test(data1)#独立性检验函数chisq.test()

 

shapiro.test(data$val)#shaprio.test()正态性检验函数

 

bartlett.test(val~trt,data)#方差齐性检验函数bartlett.test()

 

可以看到三大参数检验P>0.05,即接受原假设，数据满足独立性，正态性和方差齐性。

 第三步：单因素方差分析和多重比较

单因素方差分析

aov<-aov(val~trt,data=data)#单因素方差分析
summary(aov)

P<0.05,拒绝原假设，各组试验结果均值存在差异，即添加不同浓度的氮肥对植物根际土壤的硝态氮浓度存在影响差异。

多重比较：介绍两种函数LSD.test()和Duncan.test（）

library(agricolae)
data2<-LSD.test(aov,'trt',p.adj = 'bonferroni')#'bonferroni'对P值进行修正
print(data2$groups)
plot(data2)

 

 

 可以看到不同组别之间的显著性差异。

data3<-duncan.test(aov,'trt')#duncan.test()函数
data3$groups
plot(data3)

第四步：提取关键数据并处理（以Duncan多重比较为例）

mean<-aggregate(data$val,by=list(trt=data$trt),FUN=mean)#计算均值
sd<-aggregate(data$val,by=list(trt=data$trt),FUN=sd)#计算标准偏差
se<-sd$x/sqrt(nrow(data)/nrow(sd))                  #计算标准误差
new<-data.frame(trt=mean$trt,mean=mean$x,sd=sd$x,se) #对计算后的数据重构
label<-data.frame(mean=data3$groups$val,label=data3$groups$groups) #添加多重比较结果
total<-merge(new,label,by="mean")                                  #合并数据
total

 第五步：根据结果画图

p<-ggplot(total,aes(trt,mean))+geom_bar(stat="identity",fill="grey50",width = 0.5)
p1<-p+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),width=0.1)+geom_text(aes(y=mean+se,label=label,vjust=-0.5,hjust="center"))

p2<-p1+theme(panel.grid.major = element_blank(),panel.grid.minor=element_blank(),panel.background = element_blank(),axis.line=element_line(colour="black",size=0.6))

p3<-p2+labs(x="不同的氮肥添加浓度",y="植物根际土壤硝态氮浓度")+ylim(0,200)
p3

 end~