方差分析(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.小結
