【數據分析 R語言實戰】學習筆記 第八章 方差分析與R實現

方差分析泛應用於商業、經濟、醫學、農業等諸多領域的數量分析研究中。例如商業廣告宣傳方面，廣告效果可能會受廣告式、地區規模、播放時段、播放頻率等多個因素的影響，經過方差分析研究衆多因素中，哪些是主要的以及如何產生影響等。而在經濟管理中，方差分析經常使用於分析變量之間的關係，如人民幣匯率對股票收益率的影響、存貸款利率對債券市場的影響，等等。

協方差是在方差分析的基礎上，綜合迴歸分析的方法，研究如何調節協變量對因變量的影響效應，從而更加有效地分析實驗處理效應的一種統計技術。

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, ...)
例:
某製造商僱用了來自三所本地大學的僱員做爲管理人員。最近，公司的人事部門已經收集信息並考覈了年度工做成績。從三所大學來的僱員中隨機地抽取了三個獨立樣本，樣本量分別爲七、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值很是小，說明結果很是顯著，應該拒絕原假設，認爲各因素在不一樣水平下的試驗結果有顯著差異，即三種肥料對蘋果產量有很大的影響。