目錄
- 七、方差分析
-
- 7.1 單因子方差分析
-
- 7.1.1 數學模型
-
- 方差分析表
- 7.1.2 均值的多重比較
-
- 多重 t 檢驗方法
-
- ==差異顯著與顯著差異的區分==
- ==置信區間的易混淆==
- 7.1.3 同時置信區間：Tukey法
- 7.1.4 方差齊性檢驗
-
- Barlett檢驗
- Levene檢驗
- 注意
- 7.2 雙因子方差分析
-
- 7.2.1 無互動作用的方差分析
- 7.2.2 有互動作用的方差分析
- 7.3 協方差分析
七、方差分析
方差分析的主要工作就是将觀測資料的總變異按照變異原因的不同分解為因子效應與試驗誤差，并對其做出數量分析，比較各種原因在總變異中所占的重要程度，以此作為進一步統計推斷的依據。
前提條件：獨立、正态、方差齊性
7.1 單因子方差分析
7.1.1 數學模型
設試驗僅有一個因子（因素）A，其有r個水準A1,A2,…,Ar，現在水準Ai下進行ni次獨立觀測，得到觀測資料Xij，j=1,2,…,ni，i=1,2,…,r
SST 總離差平方和（總變差）：Xij與總平均值之差的平方和，描繪所有資料的離散程度。
SSE 誤差平方和（組内平方和）：對固定的水準i，觀測值Xij之間的差異大小的度量。
SSA 因素A的效應平方和（組間平方和）：因子A各水準下的樣本均值和總平均值的平方和。
formula 是方差分析的公式，在單因素方差分析中表示為 x~A, data 是資料框。
例： 對除雜方法進行選擇。在實驗中選用5種不同的除雜方法，每種方法做4次試驗，即重複4次，如下表所示。
除雜方法Ai 除雜量Xij 平均量 A1 25.6 22.2 28.0 29.8 26.4 A2 24.4 30.0 29.0 27.5 27.7 A3 25.0 27.7 23.0 32.2 27.0 A4 28.8 28.0 31.5 25.9 28.6 A5 20.6 21.2 22.0 21.2 21.3
> x<- c(25.6, 22.2, 28.0, 29.8,24.4 ,30.0 ,29.0, 27.5,25.0 ,27.7, 23.0 ,32.2,28.8, 28.0, 31.5 ,25.9,20.6, 21.2, 22.0 ,21.2)
> A<-factor(rep(1:5,each=4))
> miscellany<-data.frame(x,A)
> aov.mis<-aov(x~A,data=miscellany)
> summary(aov.mis)
Df Sum Sq Mean Sq F value Pr(>F)
A 4 132 33.0 4.31 0.016 *
Residuals 15 115 7.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Df 表示自由度，Sum Sq 表示平方和，Mean Sq表示均方和，F value 表示F檢驗統計量的值，Pr(>F)表示檢驗p的值， A 就是因素， Residuals 為殘差。
可以看出F=4.31>F0.05(5-1,20-5)=3.06，或p=0.0016<0.05說明有理由拒絕原假設，即認為五種除雜方法的差異顯著。
方差分析表
據上述結果可以填寫下面的方差分析表：
| 方差來源 | 自由度df | 平方和SS | 均方和MS | F比 | p值 |
|---|---|---|---|---|---|
| 因素A | 4 | 132 | 33.0 | 4.31 | 0.016 |
| 殘差E | 15 | 115 | 7.7 | ||
| 總和T | 19 | 247 |
從上圖也可以看出5種除雜方法産生的除雜量有顯著差異，特别是第5種與前4種，而方法1和3，方法2和4的差異不明顯。
7.1.2 均值的多重比較
之前在進行方差分析後發現各效應的均值之間有顯著差異，隻能知道有某些均值彼此不同，但不知道具體是哪些均值不同。
多重 t 檢驗方法
當多次重複使用該方法時，會增大犯第一類錯誤的機率，進而使得結論不一定可靠。是以在進行較多次重複比較時，要對p值進行調整：
p是p值構成的向量，method是修正方法。
當比較次數較多時，Bonferroni方法的效果較好。
得到多重比較的p值：
x是響應變量構成的向量，g是分組向量（因子）, p.adjust.methods是上面提到的調整p值的方法，p.adjust.methods=none表示不做任何調整。預設值按Holm方法進行調整。
例：對上例作均值的多重比較，進一步檢驗。
下面用三種方法進行多重比較：
1）不對p作出調整
> pairwise.t.test(x,A,p.adjust.method = 'none')
Pairwise comparisons using t tests with pooled SD
data: x and A
1 2 3 4
2 0.509 - - -
3 0.773 0.707 - -
4 0.289 0.679 0.434 -
5 0.019 0.005 0.010 0.002
P value adjustment method: none
2）按預設值”holm“對p值進行調整
> pairwise.t.test(x,A)
Pairwise comparisons using t tests with pooled SD
data: x and A
1 2 3 4
2 1.00 - - -
3 1.00 1.00 - -
4 1.00 1.00 1.00 -
5 0.13 0.04 0.08 0.02
P value adjustment method: holm
3）按”bonferroni“對p值進行調整
> pairwise.t.test(x,A,p.adjust.method = "bonferroni")
Pairwise comparisons using t tests with pooled SD
data: x and A
1 2 3 4
2 1.00 - - -
3 1.00 1.00 - -
4 1.00 1.00 1.00 -
5 0.19 0.05 0.10 0.02
P value adjustment method: bonferroni
從輸出結果可以看出，做調整後p值增大，在一定程度上克服了多重t檢驗的缺點。
差異顯著與顯著差異的區分
由于常用“顯著”來表示P值大小，是以P值最常見的誤用是把統計學上的顯著與臨床或實際中的顯著差異相混淆，即混淆“差異具有顯著性”和“具有顯著差異”二者的意思。其實，前者指的是p<=0.05，即說明有充分的理由認為比較的二者來自同一總體的可能性不足5%，因而認為二者确實有差異，下這個結論出錯的可能性<=5%。而後者的意思是二者的差别确實很大。舉例來說，4和40的差别很大，因而可以說是“有顯著差異”，而4和4.2差别不大，但如果計算得到的P值<=0.05，則認為二者“差别有顯著性”，但是不能說“有顯著差異”。
————————————————
原文連結：https://blog.csdn.net/yu1581274988/article/details/117295802
置信區間的易混淆
參數95%的置信度在區間A的意思是：
正确：采樣100次計算95%置信度的置信區間，有95次計算所得的區間包含真實值。
錯誤：采樣100次，有95次真實值落在置信區間。
真實值不會變，變的是置信區間。
————————————————
原文連結：https://blog.csdn.net/yimingsilence/article/details/78084810
7.1.3 同時置信區間：Tukey法
對效應之差αi-αj做出置信區間，由此了解哪些效應不相等。
x是方差分析的對象, which是給出需要計算比較區間的因子向量, ordered為T則因子的水準遞增排序，進而使得因子間差異均以正值出現。
例：某商店以各自的銷售方式賣出新型手表，連續四天手表的銷量如下表所示。試考察銷售方式之間是否有顯著差異？
銷售方式 銷售量資料 A1 23 19 21 13 A2 24 25 28 27 A3 20 18 19 15 A4 22 25 26 23 A5 24 23 26 27
#生成資料
> sales<-data.frame(x=c(23, 19, 21, 13,24 ,25, 28, 27,20 ,18 ,19 ,15,22, 25 ,26 ,23,24 ,23 ,26, 27),A=factor(rep(1:5,c(4,4,4,4,4))))
#進行方差分析
> summary(aov(x~A,sales))
Df Sum Sq Mean Sq F value Pr(>F)
A 4 213 53.2 7.98 0.0012 **
Residuals 15 100 6.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#可見，不同銷售方式有差異，求均值之差的同時置信區間
> TukeyHSD(aov(x~A,sales))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = x ~ A, data = sales)
$A
diff lwr upr p adj
2-1 7 1.3622 12.638 0.0120
3-1 -1 -6.6378 4.638 0.9806
4-1 5 -0.6378 10.638 0.0945
5-1 6 0.3622 11.638 0.0344
3-2 -8 -13.6378 -2.362 0.0042
4-2 -2 -7.6378 3.638 0.8062
5-2 -1 -6.6378 4.638 0.9806
4-3 6 0.3622 11.638 0.0344
5-3 7 1.3622 12.638 0.0120
5-4 1 -4.6378 6.638 0.9806
7.1.4 方差齊性檢驗
檢驗資料在不同水準下方差是否相同。
Barlett檢驗
H0：各因子水準下的方差相同，H1：各因子水準下的方差不齊。
bartlett.test(x, g, ...)
bartlett.test(formula, data, subset, na.action, ...)
x是由資料構成的向量或清單，g是由因子構成的向量（當x是清單時，此項無效）；
formula是方差分析公式，data是資料框。
Levene檢驗
levene.test(x,group)
leveneTest(y, group, center=median, ...)
x是資料構成的向量，g是由因子構成的向量。
例：對上例資料作方差齊性檢驗。
1）用 Barlett檢驗
> bartlett.test(x~A,data=sales)
Bartlett test of homogeneity of variances
data: x by A
Bartlett's K-squared = 3.7, df = 4, p-value = 0.4
由于p=0.4>0.05，故接受原假設，認為各處理組資料是等方差的。
2）用 Levene檢驗
> leveneTest(sales$x,sales$A)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.82 0.53
15
注意
- 方差分析模型可視為一種特殊的線性模型，是以可以使用線性模型函數
lm()
anova()
aov(formula)等價于anova(lm(formula))
- 單因子方差分析還可以使用函數
oneway.test()
var.equal=TRUE
var.equal=FALSE
- 當各水準下的分布未知時，則采用Kruskal-Wallis秩和檢驗進行方差分析。
7.2 雙因子方差分析
7.2.1 無互動作用的方差分析
仍可以用函數
aov()
，方差模型公式為
x~A+B
，加号表示兩個因素具有可加性。
例：原來檢驗果汁中含鉛量有三種方法A1，A2，A3，現研究出另一種快速檢驗法A4，能否用A4替代前三種方法，需要通過試驗考察。觀察的對象是果汁，不同的果汁當作不同的水準：B1為蘋果，B2為葡萄汁，B3為蕃茄汁，B4為蘋果汁，B5為橘子汁，B6為鳳梨檸檬汁。現進行雙因素交錯搭配試驗，即用四種方法同時檢驗每一種果汁，其檢驗結果如下表所示。問因素A（檢驗方法）和B（果汁品衆）對果汁的含鉛量是否有顯著影響？
A B1 B2 B3 B4 B5 B6 xi· A1 0.05 0.46 0.12 0.16 0.84 1.30 2.93 A2 0.08 0.38 0.40 0.10 0.92 1.57 3.45 A3 0.11 0.43 0.05 0.10 0.94 1.10 2.73 A4 0.11 0.44 0.08 0.03 0.93 1.15 2.74 x·j 0.35 1.71 0.65 0.39 3.63 5.12 X··=11.85
> juice<-data.frame(x=c(0.05,0.46,0.12,0.16,0.84,1.30,0.08,0.38,0.40,0.10,0.92,1.57,0.11,0.43,0.05,0.10,0.94,1.10,0.11,0.44,0.08,0.03,0.93,1.15),A=gl(4,6),B=gl(6,1,24))
> #gl(n,k,length=n*k,labels=1:n,ordered=FALSE) n是水準數，k是每一水準上的重複次數，length是總觀測值數，ordered指明各水準是否先排序。
> juice.aov<-aov(x~A+B,data=juice)
> summary(juice.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 3 0.06 0.019 1.63 0.22
B 5 4.90 0.980 83.98 2e-10 ***
Residuals 15 0.18 0.012
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> #p說明B對含鉛量有顯著影響，而沒有充分理由說明A對鉛含量有顯著影響。
> bartlett.test(x~A,data=juice)
Bartlett test of homogeneity of variances
data: x by A
Bartlett's K-squared = 0.27, df = 3, p-value = 1
> bartlett.test(x~B,data=juice)
Bartlett test of homogeneity of variances
data: x by B
Bartlett's K-squared = 17, df = 5, p-value = 0.004
> #因素B不滿足方差齊性要求。
7.2.2 有互動作用的方差分析
仍可以用函數
aov()
，方差模型公式為
x~A+B+A:B
或
x~A*B
例：有一個關于檢驗毒品強弱的試驗，給48隻老鼠注射三種毒藥（因素A），同時有4種治療方案（因素B），這樣的試驗在每一種因素組合下都重複四次測試老鼠的存活時間，資料如下表所示。試分析毒藥和治療方案以及它們的互動作用對老鼠存活時間有無顯著影響？
A B C D 1 0.31 0.45 0.46 0.43 0.82 1.10 0.88 0.72 0.43 0.45 0.63 0.76 0.45 0.71 0.66 0.62 2 0.36 0.29 0.40 0.23 0.92 0.61 0.49 1.24 0.44 0.35 0.31 0.40 0.56 1.02 0.71 0.38 3 0.22 0.21 0.18 0.23 0.30 0.37 0.38 0.29 0.23 0.25 0.24 0.22 0.30 0.36 0.31 0.33
#建立資料框
> rats<-data.frame(x=c(0.31 ,0.45, 0.46 ,0.43,0.82 ,1.10, 0.88, 0.72,0.43, 0.45, 0.63, 0.76,0.45, 0.71 ,0.66 ,0.62,0.36, 0.29 ,0.40 ,0.23,0.92, 0.61, 0.49 ,1.24,0.44 ,0.35, 0.31, 0.40,0.56, 1.02, 0.71, 0.38,0.22 ,0.21 ,0.18 ,0.23,0.30, 0.37 ,0.38, 0.29,0.23 ,0.25, 0.24 ,0.22,0.30 ,0.36, 0.31, 0.33),A=gl(3,16,48,labels = c("1","2","3")),B=gl(4,4,48,labels = c("A","B","C","D")))
> par(mfrow=c(1,2)) #修改圖的布局
> plot(x~A+B,data=rats) #顯示各因素水準均有較大差異
用函數
interaction.plot()
作出互動效應圖，以考查因素之間互動作用是否存在。
> with(rats,interaction.plot(A,B,x,trace.label = "B"))
> with(rats,interaction.plot(B,A,x,trace.label = "A"))
兩圖中的曲線并沒有明顯相交的情況出現，是以初步認為兩個因素沒有互動作用。進一步用aov()函數檢驗。
> rats.aov<-aov(x~A*B,data = rats)
> summary(rats.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 2 1.033 0.517 23.22 3.3e-07 ***
B 3 0.921 0.307 13.81 3.8e-06 ***
A:B 6 0.250 0.042 1.87 0.11
Residuals 36 0.801 0.022
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
發現毒藥、治療方案對老鼠的治療時間有顯著影響，而互動作用的影響不顯著。
最後檢驗因素A、B下的資料是否滿足方差齊性的要求：
> bartlett.test(rats$x,rats$A)
Bartlett test of homogeneity of variances
data: rats$x and rats$A
Bartlett's K-squared = 26, df = 2, p-value = 2e-06
> bartlett.test(rats$x,rats$B)
Bartlett test of homogeneity of variances
data: rats$x and rats$B
Bartlett's K-squared = 13, df = 3, p-value = 0.004
發現p值均小于0.05，不滿足要求，這與一開始的箱線圖是一緻的。
7.3 協方差分析
協方差分析ANCOVA：是将線性回歸分析與方差分析結合起來的一種統計方法。基本思想是：将一些對響應變量Y有影響的變量（指未知或難以控制的因素）看作協變量，建立響應變量Y随協變量X變化的線性回歸關系，并利用這種回歸關系把X值化為相等後再對各處理組Y的修正均值間的差别進行假設檢驗。其實質就是從Y的總平方和中扣除X對Y的回歸平方和，對殘差平方和作進一步分解後，再進行方差分析，以更好地評價這種處理的效應。
在比較兩組或多組均數間差别的同時扣除或均衡這些不可控因素的影響
ancova(formula, data.in = NULL, ...,
x, groups, transpose = FALSE,
display.plot.command = FALSE,
superpose.level.name = "superpose",
ignore.groups = FALSE, ignore.groups.name = "ignore.groups",
blocks, blocks.pch = letters[seq(levels(blocks))],
layout, between, main,
pch=trellis.par.get()$superpose.symbol$pch)
formula是協方差分析的公式，data.in是資料框，x是協方差分析中的協變量，在作圖時若formula中沒有x則需要指出，groups是因子。
例：為研究A,B,C，三種飼料對豬的催肥效果，用每種飼料喂養8頭豬一段時間，測得每頭豬的初始重量X和增重Y。試分析三種飼料對豬的催肥效果是否相同？
A B C X1 Y1 X2 Y2 X3 Y3 1 15 85 17 97 22 89 2 13 83 16 90 24 91 3 11 65 18 100 20 83 4 12 76 18 95 23 95 5 12 80 21 103 25 100 6 16 91 22 106 27 102 7 14 84 19 99 30 105 8 17 90 18 94 32 110
飼料是認為可以控制的定性因素；豬的初始重量是難以控制的定量因子，為協變量X；實驗的觀察名額是豬的增量，為響應變量Y；各組的增重由于受豬的原始體重影響，不能直接進行方差分析，需進行協方差分析。
#建立資料集
> feed<-rep(c('A','B','C'),each=8)
> X<-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)
> Y<-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)
> pig<-data.frame(feed,X,Y)
1）認為在三種不同飼料喂養下，豬的初試體重不同，但增長速度相同
> ancova(Y~X+feed,data=pig)
Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
X 1 1621 1621 142.4 1.5e-10 ***
feed 2 707 354 31.1 7.3e-07 ***
Residuals 20 228 11
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可見豬的初始體重和增長速度對豬的增重都有顯著差異。
2）認為在三種不同飼料喂養下，豬的初始體重和增長速度都不同
> ancova(Y~X*feed,data=pig)
Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
X 1 1621 1621 162.49 1.9e-10 ***
feed 2 707 354 35.44 5.7e-07 ***
X:feed 2 48 24 2.41 0.12
Residuals 18 180 10
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
從兩個圖比較初始體重對增長的檢驗發現，三種飼料對豬的催肥效果相同，豬的初始體重對其影響不大。