方差分析泛應用於商業、經濟、醫學、農業等諸多領域的數量分析研究中。例如商業廣告宣傳方面，廣告效果可能會受廣告式、地區規模、播放時段、播放頻率等多個因素的影響，通過方差分析研究眾多因素中，哪些是主要的以及如何產生影響等。而在經濟管理中，方差分析常用於分析變量之間的關系，如人民幣匯率對股票收益率的影響、存貸款利率對債券市場的影響，等等。
協方差是在方差分析的基礎上，綜合回歸分析的方法，研究如何調節協變量對因變量的影響效應，從而更加有效地分析實驗處理效應的一種統計技術。
8.1單因素方差分析及R實現
(1)正態性檢驗
對數據的正態性，利用Shapiro-Wilk正態檢驗方法(W檢驗)，它通常用於樣本容量n≤50時，檢驗樣本是否符合正態分布。
R中，函數shapiro.test()提供了W統計量和相應P值，所以可以直接使用P值作為判斷標准，其調用格式為shapiro.test(x)，參數x即所要檢驗的數據集，它是長度在35000之間的向量。
例:
某銀行規定VIP客戶的月均賬戶余額要達到100萬元，並以此作為比較各分行業績的一項指標。這里分行即因子，賬戶余額是所要檢驗的指標，先從三個分行中，分別隨機抽取7個VIP客戶的賬戶。為了用單因素方差分析判斷三個分行此項業績指標是否相同，首先對二個分行的賬戶余額分別進行正態檢驗。
> x1=c(103,101,98,110,105,100,106) > x2=c(113,107,108,116,114,110,115) > x3=c(82,92,84,86,84,90,88) > shapiro.test(x1) Shapiro-Wilk normality test data: x1 W = 0.97777, p-value =0.948 > shapiro.test(x2) Shapiro-Wilk normality test data: x2 W = 0.91887, p-value =0.4607 > shapiro.test(x3) Shapiro-Wilk normality test data: x3 W = 0.95473, p-value =0.7724
P值均大於顯著性水平a=0.05,因此不能拒絕原假設，說明數據在因子A的三個水平下都
是來自正態分布的。
(2)方差齊性檢驗
方差分析的另一個假設：方差齊性，需要檢驗不同水平卜的數據方差是否相等。R中最常用的Bartlett檢驗,bartlett.test()調用格式為
bartlett.test(x，g…)
其中，參數X是數據向量或列表(list) ; g是因子向量，如果X是列表則忽略g.當使用數據集時，也通過formula調用函數:
bartlett.test(formala, data, subset，na.action…)
formula是形如lhs一rhs的方差分析公式;data指明數據集:subset是可選項，可以用來指定觀測值的一個子集用於分析:na.action表示遇到缺失值時應當采取的行為。
續上例:
> x=c(x1,x2,x3) > account=data.frame(x,A=factor(rep(1:3,each=7))) > bartlett.test(x~A,data=account) Bartlett test of homogeneity of variances data: x by A Bartlett's K-squared = 0.13625, df = 2, p-value = 0.9341
由於P值遠遠大於顯著性水平a=0.05，因此不能拒絕原假設，我們認為不同水平下的數據是等方差的。
8.1.2單因素方差分析
R中的函數aov()用於方差分析的計算，其調用格式為:
aov(formula, data = NULL, projections =FALSE, qr = TRUE,contrasts = NULL, ...)
其中的參數formula表示方差分析的公式，在單因素方差分析中即為x~A ; data表示做方差分析的數據框:projections為邏輯值，表示是否返回預測結果:qr同樣是邏輯值，表示是否返回QR分解結果，默認為TRUE; contrasts是公式中的一些因子的對比列表。通過函數summary()可列出方差分析表的詳細結果。
上面的例子已經對數據的正態性和方差齊性做了檢驗，接F來就可以進行方差分析:
> a.aov=aov(x~A,data=account) > summary(a.aov) Df Sum Sq Mean Sq F value Pr(>F) A 2 2315 1158 82.68 8.46e-10 *** Residuals 18 252 14 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plot(account$x~account$A)
Levene檢驗
Levene檢驗，它既可以用於正態分布的數據，也可用於非正態分布的數據或分布不明的數據，具有比較穩健的特點，檢驗效果也比較理想。
R的程序包car中提供了Levene檢驗的函數levene.test()
> library(car) > levene.test(account$x,account$A) Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F) group 2 0.0426 0.9584 18
由於p值大於a=0.05，不能拒絕原假設，我們認為不同水平下的數據是等方差的。
8.1.3多重t檢驗
單因素方差分析是從總體的角度上說明各效應的均值之間存在顯著差異，但具體哪些水平下的均值存在較人差異無從得知，所以我們要對每一對樣本均值進行一一比較，即要進行均值的多重比較。
> p.adjust.methods [1] "holm" "hochberg" "hommel" "bonferroni" "BH" [6] "BY" "fdr" "none" > attach(account) > pairwise.t.test(x,A,p.adjust.method="bonferroni") Pairwise comparisons using t tests with pooled SD data: x and A 1 2 2 0.0013 - 3 3.9e-07 6.5e-10 P value adjustment method: bonferroni
經過修正后的p值比原來會增大很多，這在一定程度上克服了多重t檢驗增加犯第一類錯誤的
概率的缺點。從檢驗結果來看，樣本兩兩之問t檢驗的p值都很小，說明幾個樣本之間差異明顯。
8.1.4Kruskal-Wallis秩和檢驗
R內置函數kruskal.test()可以完成Kruskal-Wallis秩和檢驗，使用如下:
kruskal.test(x, ...)
kruskal.test(x, g, ...)
kruskal.test(formula, data, subset,na.action, ...)
例:
某制造商雇用了來自三所本地大學的雇員作為管理人員。最近，公司的人事部門已經收集信息並考核了年度工作成績。從三所大學來的雇員中隨機地抽取了三個獨立樣本，樣本量分別為7、6, 7，數據如表所示。制造商想知道來自這三所不同的大學的雇員在管理崗位上的表現是否有所不同，我們通過Kruskal-Wallis秩和檢驗來得到結論。
>data=data.frame(x=c(25,70,60,85,95,90,80,60,20,30,15,40,35,50,70,60,80,90,70,75),g=factor(rep(1:3,c(7,6,7)))) > kruskal.test(x~g, data=data) Kruskal-Wallis rank sum test data: x by g Kruskal-Wallis chi-squared = 8.9839, df = 2, p-value = 0.0112
檢驗的結果為P=0.0112<0.05，因此拒絕原假設，說明來自這三個不同的大學的雇員在管理崗位上的表現有比較顯著的差異。
8.2雙因素方差分析及R實現
8.2.1無交互作用的分析
例：
某商品在不同地區、不同包裝的銷售數據
首先為了建立數據集，引入生成因子水平的函數g1()，其調用格式為：
gl(n, k, length=n*k,labels=1:n,ordered=FALSE)
n是因子的水平個數;k表示每一水平上的重復次數;length=n*k表示總觀測數;可通過參數labels對因子的不同水平添加標簽;ordered為邏輯值，指示是否排序。
> x=c(20,12,20,10,14,22,10,20,12,6,24,14,18,18,10,16,4,8,6,18,26,22,16,20,10) > sales=data.frame(x,A=gl(5,5),B=gl(5,1,25)) > sales$B [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 12 3 4 5 Levels: 1 2 3 4 5
分析前先對因素A和B作方差齊性檢驗，使用函數bartlett.test()
> bartlett.test(x~A,data=sales) Bartlett test of homogeneity of variances data: x by A Bartlett's K-squared =0.66533, df = 4, p-value = 0.9555 > bartlett.test(x~B,data=sales) Bartlett test of homogeneity of variances data: x by B Bartlett's K-squared =1.2046, df = 4, p-value = 0.8773
因素A和B的P值都遠大於0.05的顯著性水平，不能拒絕原假設，說明因素A, B的各水平是滿足方差齊性的。這時再進行雙因素方差分析，輸入指令
> sales.aov=aov(x~A+B,data=sales) > summary(sales.aov) Df Sum Sq Mean Sq F valuePr(>F) A 4 199.4 49.84 2.303 0.1032 B 4 335.4 83.84 3.874 0.0219 * Residuals 16 346.2 21.64 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1
檢驗的結論:因素B的P值=0.0219<0.05，拒絕原假設，說明銷售地區對飲料的銷售量有顯著影響;而因素A的P值=0.1032>0.05，不能拒絕原假設，因此沒有充分的理由可以說明包裝方式對銷售有明顯影響。
8.2.2有交互作用的分析
R仍然用函數aov()作雙因素方差分析，只需將formula改為x~A+B+A:B或x~A*B的形式即可。
例：
不同路段和不同時段的行車時間數據
首先構造數據集，對因素A和B作方差齊性檢驗，利用函數bartlett.test()
> time=c(25,24,27,25,25,19,20,23,22,21,29,28,31,28,30,20,17,22,21,17,18,17,13,16,12,22,18,24,21,22)
> traffic=data.frame(time,A=gl(2,15,30),B=gl(3,5,30,labels=c("I","II","III")))
> bartlett.test(time~A,data=traffic)
Bartlett test of homogeneity of variances
data: time by A
Bartlett's K-squared =0.053302, df = 1, p-value = 0.8174
> bartlett.test(time~B,data=traffic)
Bartlett test of homogeneity of variances
data: time by B
Bartlett's K-squared =0.57757, df = 2, p-value = 0.7492
檢驗結果的P值均遠大於顯著性水平0.05，說明兩個因素下的各水平都滿足方差齊性的要求，可以進一步做方差分析。畫圖來觀察一下數據的特點，首先是箱線圖。
> op=par(mfrow=c(1,2)) #分割圖形區域 > plot(time~A+B,data=traffic) Hit <Return> tosee next plot:
從圖形上單獨觀察時段和路段對行車時間的影響，可以發現因素的不同水平還是有明顯差別的。為了考察因素間的交互作用是否存在，利用函數interaction.plot()繪制交互效應圖：
interaction.plot(x.factor, trace.factor,response, fun = mean,type = c("l","p", "b", "o", "c"), legend = TRUE,trace.label =deparse(substitute(trace.factor)),fixed = FALSE,xlab =deparse(substitute(x.factor)),ylab = ylabel,ylim = range(cells, na.rm =TRUE),lty = nc:1, col = 1, pch =c(1:9, 0, letters),xpd = NULL, leg.bg =par("bg"), leg.bty = "n",
xtick = FALSE, xaxt = par("xaxt"),axes = TRUE,...)
x.factor表示橫軸的因子
trace.factor表示分類繪圖的因子
response是數值向量，要輸入響應變量
fun表示匯總數據的方式，默認為計算每個因子水平下的均值
type指定圖形類型
legend是邏輯值，指示是否生成圖例
trace.label給出圖例中的標簽。
> attach(traffic) > interaction.plot(A,B,time,legend=F) > interaction.plot(B,A,time,legend=F)
曲線均沒有相交，所以可以初步判斷兩個因素之間應該沒有交互作用。用方差分析進行確認：
> traf.aov=aov(time~A*B,data=traffic) > summary(traf.aov) Df Sum Sq Mean Sq F value Pr(>F) A 1 313.63 313.63 84.766 2.41e-09 *** B 2 261.60 130.80 35.351 7.02e-08 *** A:B 2 6.67 3.33 0.901 0.42 Residuals 24 88.80 3.70 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根據檢驗結果的P值作判斷:引素A時段和B路段對行車時間有顯著影響;而交互作用A:B的P值=0.42>0.05 ,因此不能拒絕原假設H0，說明兩個因素間沒有明顯的交互效應。
8.3協方差分析及R實現
為了提高試驗的精確性和准確性，我們對除研究因素以外的一切條件都需要采取有效措施嚴加控制，使它們在因素的不同水平間盡量保持一致，這叫做試驗控制。但當我們進行試驗設計時，即使做出很大努力控制，也經常會碰到試驗個體的初始條件不同的情況，如果不考慮這些因素有可能導致結果失真。如果考慮這些不可控的因素，這種方差分析就叫做協方差分析，其是將回歸分析和方差分析結合在一起的方法。它的基本原理如下:將一些對響應變量Y有影響的變量X(未知或難以控制的因素)看作協變量，建立響應變量Y隨X變化的線性回歸分析，從Y的總的平方和中扣除X對Y的回歸平方和，對殘差平方和作進一步分解后再進行方差分析。
例：
施用3種肥料的蘋果產量
> Weight_Initial=c(15,13,11,12,12,16,14,17,17,16,18,18,21,22,19,18,22,24,20,23,25,27,30,32) > Weight_Increment=c(85,83,65,76,80,91,84,90,97,90,100,95,103,106,99,94,89,91,83,95,100,102,105,110) > feed=gl(3,8,24) > data_feed=data.frame(Weight_Initial,Weight_Increment,feed) > library(HH) > m=ancova(Weight_Increment~Weight_Initial+feed,data=data_feed) > summary(m) Df Sum Sq Mean Sq F value Pr(>F) Weight_Initial 1 1621.1 1621.1 142.44 1.50e-10 feed 2 707.2 353.6 31.07 7.32e-07 Residuals 20 227.6 11.4 Weight_Initial *** feed *** Residuals --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
協方差分析的P值非常小，說明結果非常顯著，應該拒絕原假設，認為各因素在不同水平下的試驗結果有顯著差別，即三種肥料對蘋果產量有很大的影響。
