R的方差检验

来源:互联网 时间:2017-09-19

方差分析(Analysis of Variance,简称ANOVA),又称“ 变异数分析”,是R.A.Fisher发明的,用于两个及两个以上 样本均数差别的 显著性检验。 由于各种因素的影响,研究所得的数据呈现波动状。造成波动的原因可分成两类,一是不可控的随机因素,另一是研究中施加的对结果形成影响的可控因素。方差分析是从观测变量的方差入手,研究诸多 控制变量中哪些变量是对观测变量有显著影响的变量。

方差分析有两个重要的前提:

  • 控制变量不同水平下观测变量的总体分布为正态分布
  • 控制变量不同水平下观测变量的总体具有相同的方差

1.单因素方差分析

1.1概述

一个控制变量的不同水平是否对观测变量产生了显著影响。

F>>1--控制变量对观测值产生显著影响。

F分布密度曲线

x<-rf(1000,df1[i],df2[i])

set.seed(12345)

x<-rnorm(1000,0,1)

Ord<-order(x,decreasing=FALSE)

#order()的返回值是对应“排名”的元素所在向量中的位置

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)

 

1.2数学模型

1.3 R程序

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)

> aov(MPG~ModelYear,data=CarData)
Call:
aov(formula = MPG ~ ModelYear, data = CarData)

Terms:
ModelYear Residuals
Sum of Squares 10401.78 13850.79
Deg. of Freedom 12 385

Residual standard error: 5.998007
Estimated effects may be unbalanced


> OneWay<-aov(MPG~ModelYear,data=CarData)
> anova(OneWay)
Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
ModelYear 12 10402 866.82 24.094 < 2.2e-16 ***
Residuals 385 13851 35.98
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


> summary(OneWay)
Df Sum Sq Mean Sq F value Pr(>F)
ModelYear 12 10402 866.8 24.09 <2e-16 ***
Residuals 385 13851 36.0
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

F=24.09,拒绝原假设,不同年份车型的MPG总体均值存在显著差异。

1.4各总体均值的可视化 

plotmeans()

install.packages("gplots")

library("gplots")

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

1.5单因素方差分析的前提假设

  • 控制变量不同水平下观测变量的总体分布为正态分布
  • 控制变量不同水平下观测变量的总体具有相同的方差

总体正态检验

Q-Qplot / K-S test

 unique(CarData$ModelYear)

###########检验方差分析的前提假设(正态性检验一)

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)

 

K-S test

总体方差齐性检验

ks.test(数值型向量名,"pnorm")

###########检验方差分析的前提假设(正态性检验二)

for(i in unique(CarData$ModelYear)){

T<-subset(CarData,CarData$ModelYear==i)

R<-ks.test(T$MPG,"pnorm")

print(R)

}

One-sample Kolmogorov-Smirnov test

data: T$MPG
D = 1, p-value < 2.2e-16
alternative hypothesis: two-sided

于正态分布有显著差异。

各方差齐性检验

#########检验方差分析的前提假设(方差齐性性检验)

library("car")

leveneTest(CarData$MPG,CarData$ModelYear, center=mean)

Levene's Test for Homogeneity of Variance (center = mean)
Df F value Pr(>F)
group 12 1.7173 0.06103 .
385
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

不能拒绝原假设,具有相同总体方差。

1.6单因素方差检验中的多重比较检验

LSD检验

 

Tukey HSD检验

对不同车型MPG的多重比较检验

##############多重比较检验

CarData<-read.table(file="CarData.txt",header=TRUE)

CarData$ModelYear<-as.factor(CarData$ModelYear)

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)

> TukeyHSD(OneWay,ordered=FALSE,conf.level=0.95)
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = MPG ~ ModelYear, data = CarData)

$ModelYear
diff lwr upr p adj
71-70 3.5603448 -1.737584676 8.8582743 0.5621668
72-70 1.0246305 -4.273298962 6.3225600 0.9999872
73-70 -0.5896552 -5.466537903 4.2872276 0.9999999
74-70 5.0140485 -0.333563533 10.3616606 0.0912557
75-70 2.5770115 -2.630295013 7.7843180 0.9127501
76-70 3.8838742 -1.170630163 8.9383786 0.3375465
77-70 5.6853448 0.387415324 10.9832743 0.0229529
78-70 6.3714559 1.382000010 11.3609119 0.0018064
79-70 7.4034483 2.152197475 12.6546991 0.0002657
80-70 16.0068966 10.755645751 21.2581474 0.0000000
81-70 12.6448276 7.393576785 17.8960784 0.0000000
82-70 14.0200222 8.854163329 19.1858812 0.0000000

71,72,73,75.....-70无差异

 

 

1.7功效分析

###################单因素方差分析的功效分析

library("pwr")

pwr.anova.test(k=13,f=0.25,sig.level=0.05,power=0.8)

Balanced one-way analysis of variance power calculation

k = 13
n = 22.15691
f = 0.25
sig.level = 0.05
power = 0.8

NOTE: n is number in each group

 

效应量是指由于因素引起的差别,是衡量处理效应大小的指标。与显著性检验不同,这些指标不受样本容量影响。它表示不同处理下的总体均值之间差异的大小,可以在不同研究之间进行比较。

#############效应量和样本量的关系曲线

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)

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)#两者结果一致

2.单因素协方差分析

2.1概述

除了控制变量外,其他变量也会对观测值产生影响。为了更准确的研究控制变量对观测量的影响,应排除其他变量的影响。协方差分析将其他变量作为协变量,并在排除协变量对观测值影响的条件下研究控制变量的影响。

2.2数学模型

2.3R函数和示例

 

 

比如,在排除weight这个协变量的影响下,检验车型对MPG的影响

################单因素协方差分析

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))

Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
weight 1 16777.8 16777.8 1665.526 < 2.2e-16 ***
ModelYear 12 3606.6 300.5 29.835 < 2.2e-16 ***
Residuals 384 3868.2 10.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Analysis of Variance Table

Response: e
Df Sum Sq Mean Sq F value Pr(>F)
CarData$ModelYear 12 0.0 0.000 0 1
Residuals 385 3868.2 10.047

多重比较检验

install.packages("effects")

library("effects")

effect("ModelYear",Result)#调整后的MPG值

plot(effect("ModelYear",Result))

tapply(CarData$MPG,INDEX=CarData$ModelYear,FUN=mean)#调整前的MPG值

前提检验

对控制变量在不同水平下与协变量的关系是否一致且无明显差异

coplot(MPG~weight|ModelYear,data=CarData)

 

3.多因素方差分析

3.1概述

3.3R函数和示例

  • x~A+B+A:B 双因素方差,其中X~A+B中A和B是不同因素的水平因子(不考虑交互作用),A:B代表交互作用生成的因子
  • x~A*B*C==A+B+C+A:B+A:C+B:C+A:B:C
  • y~(A+B+C)^2==A+B+C+A:B+A:C+B:C
  • y~. 分析除了y以外其他全部因素对观测变量的影响

注意:~右边控制变量的排列顺序很关键。 ep, A+B与B+A是不一样的

 

################多因素方差分析

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)


Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
cylinders 4 15454.8 3863.7 280.3702 <2e-16 ***
ModelYear 12 3423.7 285.3 20.7036 <2e-16 ***
cylinders:ModelYear 26 482.0 18.5 1.3451 0.1236
Residuals 355 4892.1 13.8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

交互项对其无影响--重新构造模型

############多因素方差分析的非饱和模型

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)

Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
cylinders 4 15454.8 3863.7 273.919 < 2.2e-16 ***
ModelYear 12 3423.7 285.3 20.227 < 2.2e-16 ***
Residuals 381 5374.1 14.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Analysis of Variance Table

Response: MPG
Df R Sum Sq R Mean Sq Iter Pr(Prob)
cylinders 4 8476.7 2119.17 5000 < 2.2e-16 ***
ModelYear 12 3423.7 285.31 5000 < 2.2e-16 ***
Residuals 381 5374.1 14.11
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.小结

 

相关阅读:
Top