R软件学习笔记-7（方差分析）

方差分析的基本步骤：

1、建立检验假设；

H0：多个样本总体均值相等；

H1：多个样本总体均值不相等或不全等。

检验水准为0.05。

2、 计算检验统计量F值；

3、 确定P值并作出推断结果。

基本假设：

1. 方差分析的假定条件为：

（1）各处理条件下的样本是 随机的。

（2）各处理条件下的样本是 相互独立的，否则可能出现无法解析的输出结果。

（3）各处理条件下的样本分别来自 正态分布总体，否则使用 非参数分析。

（4）各处理条件下的 样本方差相同，即具有齐效性。

应用条件：

各样本是相互独立的随机样本

1. 各样本是 相互独立的随机样本
2. 各样本均来自 正态分布总体

3. 各样本的总体方差相等，即具有 方差齐性

4.在 不满足正态性时可以用非参数检验

2. 方差分析的 假设检验

假设有K个样本，如果原 假设H0样本均数都相同， K个样本有共同的方差σ ，则 K个样本来自具有共同方差σ和相同均值的总体。

如果经过计算，组间 均方远远大于组内均方，则推翻原假设，说明样本来自不同的正态总体，说明处理造成均值的差异有统计意义。否则承认原假设，样本来自相同总体，处理间无差异。

###############F分布

set.seed(12345)
x<-rnorm(1000,0,1)
Ord<-order(x,decreasing=FALSE)
x<-x[Ord]
y<-dnorm(x,0,1)

plot(x,y,xlim=c(-1,5),ylim=c(0,2),type="l",ylab="密度",main="标准正态分布与不同自由度下的F分布密度函数",lwd=1.5)

#######不同自由度的F分布
df1<-c(10,15,30,100)
df2<-c(10,20,25,110)
for(i in 1:4){
 x<-rf(1000,df1[i],df2[i])
 Ord<-order(x,decreasing=FALSE)
 x<-x[Ord]
 y<-df(x,df1[i],df2[i])
 lines(x,y,lty=i+1)
}
legend("topright",title="自由度",c("标准正态分布",paste(df1,df2,sep="-")),lty=1:5)          ##legend函数参见后面的的文章，paste函数表示的是将两个向量粘合在一起，中间用seq做分割


单因素方差分析
##################单因素方差分析
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
aov(MPG~ModelYear,data=CarData)
OneWay<- aov(MPG~ModelYear,data=CarData)
anova(OneWay)

summary(OneWay)

均值变化折线图
##############均值变化折线图（带置信区间）
install.packages("gplots")
library("gplots")
plotmeans(MPG~ModelYear,data=CarData,p=0.95,use.t=TRUE,xlab="年代车型",ylab="平均MPG",main="不同年代车型MPG总体均值变化折线图(95%置信区间)")

正态性检验一

###########检验方差分析的前提假设（正态性检验一）
par(mfrow=c(3,5),mar=c(4,4,4,4))
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 qqnorm(T$MPG,main=paste(i,"年车型mpg Q-Q图"),cex=0.7)
 qqline(T$MPG,distribution = qnorm)
}

############或者
library("lattice")
qqmath(~MPG|ModelYear,data=CarData)

正态性检验二

###########检验方差分析的前提假设（正态性检验二）
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 R<-ks.test(T$MPG,"pnorm")
 print(R)

}

方差齐性性检验

#########检验方差分析的前提假设（方差齐性性检验）
install.packages("car")
library("car")
leveneTest(CarData$MPG,CarData$ModelYear, center=mean)

多重比较检验

##############多重比较检验
OneWay<-aov(MPG~ModelYear,data=CarData)
OneWay$coefficients
TukeyHSD(OneWay,ordered=FALSE,conf.level=0.95)
Result<-TukeyHSD(OneWay,ordered=TRUE,conf.level=0.95)
LineCol<-vector()
LineCol[Result[[1]][,4]<0.05]<-2
LineCol[Result[[1]][,4]>=0.05]<-1
par(las=2)
par(mar=c(5,8,4,2))
plot(Result,cex.axis=0.5,col=LineCol)

功效分析

###################单因素方差分析的功效分析
library("pwr")
pwr.anova.test(k=13,f=0.25,sig.level=0.05,power=0.8)

#############效应量和样本量的关系曲线
library("pwr")
ES<-seq(from=0.1,to=0.8,by=0.01)
SampleSize<-matrix(nrow=length(ES),ncol=8)
for(i in 3:10){
 for(j in 1:length(ES)){
  result<-pwr.anova.test(k=i,f=ES[j],sig.level=0.05,power=0.8)
  SampleSize[j,i-2]<-ceiling(result$n)
  }
 }
plot(SampleSize[,1],ES,type="l",ylab="效应量",xlab="样本量（每个水平）",main="单因素方差分析(Alpha=0.05,Power=0.8)")
for(i in 2:8){
 lines(SampleSize[,i],ES,type="l",col=i)
}
legend("topright",title="水平数",paste("k",3:10,sep="="),lty=1,col=1:8)

###################单因素方差分析的置换检验
install.packages("lmPerm")
library("lmPerm")
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
OneWay<-aov(MPG~ModelYear,data=CarData)
anova(OneWay)
Fit<-aovp(MPG~ModelYear,data=CarData,perm="Prob")
anova(Fit)

################单因素协方差分析
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
Result<-aov(MPG~weight+ModelYear,data=CarData)
anova(Result)

e<-CarData$MPG-Result$fitted.values   #剔除协变量影响后的残差
anova(aov(e~CarData$ModelYear))

var(CarData$MPG)*(398-1)-16777.8-3868.2

install.packages("effects")
library("effects")
effect("ModelYear",Result)
plot(effect("ModelYear",Result))
tapply(CarData$MPG,INDEX=CarData$ModelYear,FUN=mean)

#################单因素协方差分析的前提检验
coplot(MPG~weight|ModelYear,data=CarData)

单因素协方差分析的置换检验
################单因素协方差分析的置换检验
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
Result<-aov(MPG~weight+ModelYear,data=CarData)
anova(Result)
library("lmPerm")
Fit<-aovp(MPG~weight+ModelYear,data=CarData,perm="Prob")
anova(Fit)

多因素方差分析
################多因素方差分析
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
CarData$cylinders<-as.factor(CarData$cylinders)
table(CarData$cylinders)
Result<-aov(MPG~cylinders+ModelYear+cylinders:ModelYear,data=CarData)
anova(Result)

############多因素方差分析的非饱和模型
Result<-aov(MPG~cylinders+ModelYear,data=CarData)
anova(Result)

######可视化交互效应
interaction.plot(CarData$ModelYear,CarData$cylinders,CarData$MPG,type="b",main="气缸数和车型对MPG的交互效应",xlab="车型",ylab="MPG均值")

###############多因素方差分析的置换检验
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
CarData$cylinders<-as.factor(CarData$cylinders)
Result<-aov(MPG~cylinders+ModelYear,data=CarData)
anova(Result)
library("lmPerm")
Fit<-aovp(MPG~cylinders+ModelYear,data=CarData)
anova(Fit)