第四章 方差分析
方差分析是一种特殊的假设检验,是判断多组数据之间平均数差异是否显著的。对多组数据若仍用前一章中的t检验一对对比较,会大大增加犯第一类错误的概率。例如有5组数据要比较,则共需比次。若H0正确,每次接受的概率为1-α=0.95,10次都接受为0.9510≈0.60,因此α′=1-0.60=0.40,即全部比较中至少犯一次第一类错误的概率为0.40,这显然是不能接受的。方差分析则是把所有这些组数据放在一起,一次比较就对所有各组间是否有差异作出判断。如果没有显著差异,则认为它们都是相同的;如发现有差异,再进一步比较是哪组数据与其他数据不同。这样,就避免了使α大大增加的弊病。下面我们先介绍一些方差分析中要用到的术语。
1. 因素 可能影响试验结果,且在试验中被考查的原因或原因组合。有时也可称为因子。例如温度、湿度、药物种类等。
2. 水平 因素在试验或观测中所处的状态。例如温度的不同值,药物的不同浓度等。
3. 主效应 反映一个因素各水平的平均响应之差异的一种度量。一个因子第i水平上所有数据的平均与全部数据的平均之差,称为该因子第i水平的主效应
4. 交互效应 由两个或更多因素之间水平搭配而产生的差异的一种度量。
5. 处理 实验中实施的因子水平的一个组合。
6. 固定因素 该因素的水平可准确控制,且水平固定后,其效应也固定。例如温度,化学药物的浓度,动植物的品系,等等。
7. 随机因素 该因素的水平不能严格控制,或虽水平能控制,但其效应仍为随机变量。例如动物的窝别(遗传因素的组合),农家肥的效果,等等。
8. 误差 除了实验中所考虑的因素之外,其他原因所引起的实验结果的变化。它可分为系统误差和随机误差:
系统误差:误差的组成部分,在对同一被测量的多次测试中,它保持不变或按某种规律变化。它的原因可为已知,也可为未知,但均应尽量消除。
随机误差:误差的组成部分,在对同一被测量的多次测试中,它受偶然因素的影响而以不可预知的方式变化。它无法消除或修正。
§4.1 单因素方差分析
单因素方差分析是指我们需要研究的因素只有一个,这一因素可以有几个不同的水平,我们的目标就是要看看这些水平的影响是否相同。为了在有随机误差的情况下进行比较,各水平都应有一定数量的重复。
为方便表述,我们对数据给出一种固定的表示法:
a:因素的水平数
n:每一水平的重复数
xij:第i水平的第j次观察值。1≤i≤a, 1≤j≤n
, 第i水平所有观察值的和
,第i水平均值
全部观察值的和
,总平均值
,第i水平上的子样方差。
方差分析中,我们用以下的线性统计模型描述每一观察值:
xij=(+(i+(ij, i=1, 2 …… a, j=1, 2, …… n (4.1)
其中(:总平均数;(i:i水平主效应;(ij:随机误差。要求(ij ~N(0,σ2),且互相独立。注意这里要求各水平有共同的方差σ2。
单因素方差分析的目的就是检验各(i是否均相同。由于因素可分为固定因素和随机因素,它们会对方差分析的过程产生不同的影响,我们分别加以讨论。
固定因素模型:
例4.1 用4种不同的配合饲料饲养30日令的小鸡,10天后计算平均日增重,得以下数据:
表4.1 不同饲料日增重值
饲料
日增重值Xij
1
55 49 62 45 51
2
61 58 52 68 70
3
71 65 56 73 59
4
85 90 76 78 69
4种饲料的效果是否相同?
例4.1 是固定因素模型,因为在配合饲料中,每种饲料的营养成份是固定的,它的效果也应是固定的。反映到线性模型中,就是(i是常数,且可要求
, (4.2)
这种对(i的限制并没有失去一般性,这是因为根据(4.1)式,如果各(i之和H不为0,则我们可把其和数移到总平均数(中去,即令(i ′= (i –H,从而使新的(i ′之和为0。同时,也只有新的(i ′才符合前述主效应的定义。
固定模型的统计假设为:H0:αi = 0, i = 1, 2 …… a
HA: αi ≠ 0, 至少对某一i
方差分析的基本思想,就是将总变差分解为各构成部分之和,然后对它们作统计检验。即:
由于
∴ (4.3)
用符号表示,上式可写成:
SST = SSA + SSe (4.4)
其中符号的意义为:
SST:总平方和;
SSA处理间平方和;
SSe:误差平方和,或处理内平方和。
它们的自由度分别为an–1, a–1和a(n–1),即自由度也作了相应分解:
an – 1 = a –1 + a(n – 1)
令,称为误差均方;,称为处理间均方;则它们的数学期望分别为:
∵E((ij)=0, (i为常数,且
∴ 原式 =
从这两个数学期望来看,我们给MSe和MSA起的名字是有道理的。MSe的期望是σ2,即随机误差ε的方差,说明它就是随机误差的一个估计量;而MSA的期望是,除了有代表随机误差的σ2外,还有一项是各水平主效应的平方和,即它代表了各处理间差异的大小。
若H0成立,则有:αi=0,i=1,2,…… a;此时E(MSA)=σ2;若H0不成立,则E(MSA)>σ2,令
(4.5)
则当H0成立时,F ~ F(a-1, na-a);否则F值有偏大的趋势。因此可用F分布表对H0是否成立进行上单尾检验。
方差分析的计算是比较繁杂的,因此常使用计算机进行计算。公式为:
(4.6)
(4.7)
现在的计算器常有统计功能,利用这样的计算器也可大大简化计算。步骤为:
1°把每一水平视为一个小样本,先求出它们的样本均值和样本方差,即。
2°把所有视为一个样本,求出它的样本方差,则
(4.8)
3°,或 , (4.9)
现在我们来计算例4.1(使用带统计功能的计算器):
例4.1 解:用计算器求出各处理的平均数和子样方差及平均数的子样方差:
饲料
1
2
3
4
52.4
61.8
64.8
79.6
127.24
41.8
54.2
54.2
66.3
216.5
代入(4.8)、(4.9)式,得:MSA = 5×127.24 = 636.2, MSe = 216.5/4 = 54.125,
查F分布表,得:F0.95(3, 16) = 3.24, F0.99(3, 16) = 5.29
∵ F>F0.99,∴拒绝H0,差异极显著。即:这4种饲料的增重效果差异极显著。
这就是方差分析中最简单的单因素固定模型的分析方法。对固定模型来说,如果结果是差异显著,一般还应进行多重比较,具体方法稍后介绍。从这一分析过程中可以很清楚地看到方差分析的基本思想,那就是不再对数据进行一对对的比较,而是对总体的方差进行分解,首先分离出随机误差所导致的变差,然后再将处理所引起的变差与它相比较,如果处理的变差明显大于随机误差,则说明各水平间的差异不能用随机误差解释,应认为各水平间有明显差异;否则则说明各水平间的不同可以认为是随机误差引起,即各水平间没有差异。这样就对多组实验之间的差异一次完成了检验,从而避免了多次检验引起的犯错误可能大大升高的问题。下面我们再来看看如果因素的效果是随机的,对方差分析的过程将产生什么影响。
随机因素模型
例4.2 随机选取4窝动物,每窝均有4只幼仔,其出生重见表4.2。不同窝出生重差异是否显著?
表4.2 动物出生重(g)
窝别
出生重Xij
1
34.7 33.3 26.2 31.6
2
33.2 26.0 28.6 32.3
3
27.1 23.3 27.8 26.7
4
32.9 31.4 25.7 28.0
例4.2是随机因素模型,因为动物的窝别是无法控制的,也无法重复,它的效果是无法预料的。随机因素的影响首先体现在线性统计模型中,它的表达式仍为:
xij = (+(i+(ij, i=1, 2, ……a, j=1, 2, ……n
但由于各水平的效应无法预料,现在(i不再能视为常数,而是随机变量了。即:
(NID意为独立正态分布)
此时一般Σ(i=0不再成立,统计假设相应变为:
H0:=0
HA:>0
这样,当H0成立时,自然有(i =0,i=1, 2, ……a;若不成立,则作为从N()中抽取的样本,各(i不可能都相同,当然也不可能均为0。此时它们的和一般也不会是0。
对于随机模型,总平方和与自由度的分解与固定模型是相同的,因为在证明平方和分解的过程中没有用到线性统计模型,因此因素类型的变化不会影响总平方和的分解。MSe的期望也没有变,因为这些推导过程中也没有使用(i的性质。但MSA的期望变了,因为(i不再是常数,也不再为0。
由于各(i与各εij 相互独立,上式的交叉项期望为零。因此有:
原式
从上述均方期望可看出,若H0成立,仍有:
而当HA成立时,F值仍有偏大的趋势。因此仍可用F分布表作上单尾检验。但这时对结果的解释却不同了。在固定模型中,结论只适用于检查的那几个水平。而在随机模型中由于是=0,因此结论可推广到这一因素的一切水平。
现在来计算例4.2:
例4.2 解:计算各处理平均数和方差,以及平均数的方差,填入下表:
窝别
1
2
3
4
31.45
30.025
26.225
29.50
4.88
13.86
11.16
4.01
10.62
39.65
代入(4.8), (4.9)式,得
查F分布表,得:F0.95(3, 12)=3.490 ∵ F<F0.95,∴接受H0,可认为出生重无显著差异。
从上述分析过程可知,当因素从固定变为随机后,其影响主要表现在改变了统计模型中参数(i的性质,使它从常数变成了随机变量。这样一来,所有涉及(i的地方都有了明显改变,包括统计假设H0和HA,均方期望E(MSA),以及最后的解释。对单因素方差分析来说,因素类型的变化没有影响统计量的计算与检验过程,这是与两个及更多因素方差分析不同之处。另外,由于随机因素的水平不能重复,因此多重比较也就变成没有意义的了。
不等重复时的情况。
方差分析的数据都是按照精心设计的实验方案收集来的,一般来说各水平应有相同的重复数。但若实验过程中由于某种原因丢失了一个或几个数据,又无法重做实验弥补,此时就变成各水平有不同的重复数了。在这种情况下上述方差分析的方法仍然可用,但计算公式及自由度都要作相应变化。令,则总自由度变为N-1,SSA的自由度仍为 a-1, SSe的自由度变为N-a。(4.6),(4.7)式相应变为
(4.10)
(4.11)
用计算器的计算方法也应改为:
1°计算每一处理的样本方差
2°全部样本放在一起,计算总样本方差S2
3°SST =(N-1)·S2
4° (4.12)
5°SSA = SST - SSe (4.13)
多重比较
固定模型拒绝H0时,并不意味着所有处理间均存在差异。为弄清哪些处理间有差异,需对所有水平作一对一的比较,即多重比较。常用的多重比较方法有以下几种:
最小显著差数(LSD)法:实际就是用t检验对所有平均数作一对一对的检验。一般情况下各水平重复数n相等,用MSe作为的估计量,可得:
统计量为:
因此当
(4.14)
时,差异显著。t分位数的自由度df = a(n-1)。
即为最小显著差数,记为LSD。所有比较仅需计算一个LSD,应用很方便。但由于又回到了多次重复使用t检验的方法,会大大增加犯第一类错误的概率。为了克服这一缺点,人们提出了多重范围检验的思想:即把平均数按大小排列后,对离得远的平均数采用较大的临界值R。这一类的方法主要有Duncan法和Newman-Keul法。后者又称为q法。现介绍如下:
2. Duncan法: Duncan法步骤如下:
1°把需比较的a个平均数从大到小排好:
2°求出各对差值,并列成表:
表4.3 a个均值间的差值表
a
a-1
…………
3
2
1
2
…
a-2
a-1
…
…
…………
……………
3°求临界值,K=2,3……a。 (4.15)
其中α=0.05或0.01,k表示两平均数在位次上的差别,即若差为,则k=j-i+1。因此相邻二平均数k值为2,隔一个为3,余类推。,df为MSe的自由度。的值需查专门表格。最后把求得的临界列成下表:
表4.4 多重检验临界值表
K
γ0.05
R0.05
γ0.01
R0.01
2
3
…
a
γ0.05(2,df)
γ0.05(3,df)
…
γ0.05(a,df)
R2,0.05
R3,0.05
…
Ra,0.05
γ0.01(2,df)
γ0.01(3,df)
…
γ0.01(a,df)
R2,0.01
R3,0.01
…
Ra,0.01
4°对差值表采用适当的R进行比较。差值表中每条对角线上的k值是相同的,可使用同一个临界值R。差值大于 R0.05,标以“*”; 大于 R0.01则标“**”。若比较的两个水平重复数不等,设为ni, nj,则可用它们的调和平均值nij代替n。即:
,
此时
(4.16)
3.Newman-Q法。又称多重范围q检验。它的检验方法与Duncan法完全相同,只是要查不同的系数表。它的系数表称为q值表。
三种方法的比较:
比较Duncan的r值表与q值表,可知当k=2时,,此时三种检验法是相同的。当k≥3时,三种方法临界值不同,其中LSD最小,Duncan法次之,Newman-Q法最大。因此LSD法犯第一类错误概率最大,Duncan法次之, Newman-Keul法最小,可按照犯两类错误危害性大小选择适当的方法。一般来说,Duncan法最常用;若各水平均值只需与对照比较,由于比较次数较少,可考虑选用LSD法。另外,只有F检验确认各平均数间有显著差异后才可进行LSD法检验,而另两种方法则不一定,有时它们的结果也可能与F检验不一致。
例4.3 对例4.1进行多重比较。
解:前已算出:,,,,
MSe = 54.125, df = 4×(5-1)=16
最小显著差数法:
查表,得t0.975(16)=2.1199, t0.995(16)=2.9208
∴ LSD0.05 = t0.975(16)=2.1199×
= 2.1199×4.6530 = 9.8639
LSD0.01 = 2.9208×4.6530 = 13.5905
列出各水平均值的差值表:(均值已从小到大排列,不必再排)
4
3
2
1
2
3
27.2**
17.8**
14.8**
12.4*
3.0
9.4
将各差值分别与LSD0.05和LSD0.01比较,大于LSD0.05的标“*”,大于LSD0.01的标“**”。得:与其他三个均值均达差异极显著,与差异显著。
(2)Duncan法:
df=16
利用公式求各临界值:
表4.5 Duncan多重检验临界值表
K
r0.05(k, 16)
R0.05
r0.01(k, 16)
R0.01
2
3
4
3.00
3.15
3.23
9.87
10.36
10.63
4.13
4.34
4.45
13.59
14.28
14.64
列出差值表,并与临界值表中的数值进行比较:
4
3
2
1
2
3
27.2**
17.8**
14.8**
12.4*
3.0
9.4
最长的对角线上应使用k=2的临界值, 因此首先与α=0.05的临界值9.87比较, 大于9.87的则标一个“*”号;再与α=0.01的临界值13.59比较,大于13.59则再加一个“*”号。次长对角线应使用k=3的临界值,因此应先后与10.36,14.28比较,大于前者加一个“*”,大于后者再加一个“*”。第三条对角线上只有一个数27.2,它应与k=4的临界值,即10.63和14.64比较,显然它比这两个临界值都大,因此也应标上两个“*”号。这样就完成了多重比较。把这一差值表与前边最小显著差数法的差值表进行比较,可以看到它们的结果是相同的。但若比较一下两种方法的临界值,就可以发现Duncan法k=2的临界值就是最小显著差数法的临界值,而k>2的Duncan法临界值变大,但对本题来说,这种变大尚不足以改变最终的结果。
(3)Newman-Q法:
仍有: df=16 。
利用公式求各临界值:
表4.6 Newman-Q法临界值表
K
q0.05(k, 16)
Q0.05
q0.01(k, 16)
Q0.01
2
3
4
3.00
3.65
4.05
9.87
12.01
13.32
4.13
4.79
5.19
13.59
15.76
17.08
列出差值表,并与相应临界值比较:
4
3
2
1
2
3
27.2**
17.8**
14.8**
12.4*
3.0
9.4
与Duncan法同样,最长的对角线使用k=2的两个临界值,即9.87和13.59比较,大于前者加“*”,大于后者再加一个“*”;右上次长对角线用k=3,即临界值12.01和15.76;最后一条用k=4,即13.32和17.08比较。最终结果与前两种方法仍相同,但与的差12.4已接近临界值12.01。
比较三种方法,当k=2时临界值均相同,当k>2时临界值依次增大;但对本例题来说,这种增大还不足以影响最终结果。
§4.2 多因素方差分析
上一节我们讨论了最简单的方差分析——单因素方差分析的原理与方法。在实际工作中,问题常常比较复杂,要求我们同时考虑两种甚至更多因素,以及这些因素共同作用的影响。此时单因素方差分析就无能为力了,需采用两因素或更多因素方差分析。进行多因素方差分析从理论上说并无任何困难,但随着因素数的增加,普通方差分析的复杂性迅速增加,这种复杂性不仅表现在分析计算的繁复,更表现在所需实验次数呈现出几何级数的增加上。这样一来,当因素数增加到三个或三个以上时,其工作量之大常常是令人望而生畏。因此三或三因素以上方差分析较少用到;当确实需要考虑这样多因素时,我们常常转而采用一些特殊的方差分析方法,例如正交实验设计方法,有关内容我们将在第九章中介绍。由于以上原因,本节内容将主要集中在讨论两因素方差分析上。
模型类型及交互作用概念。
与单因素方差分析相比,交互作用是多因素方差分析中新的概念之一。当一个因素的效应明显地依赖于其他因素的水平时,我们称这些因素间有交互效应。例如,由于人的体质不同,药物的疗效也可能会有不同;不同的地施用同样的肥料,增产效果也有不同,等等。交互效应的有无可用一些直观方法粗略估计,例如可用图形来估计:
B1
B2
B3
A1 A2 A3
B1
B2
B3
A1 A2 A3
(a) 无交互效应 (b) 有交互效应
图4.1 交互效应示意图
图中每条曲线代表B因素的一个水平。若各曲线平行或近似平行,可认为无交互效应,否则为有交互效应。以上只是一种直观的判断,在多因素方差分析的过程中,我们对交互作用的有无也可进行统计检验。具体原理与方法我们将在下文中详细介绍。
多因素方差分析可按照不同标准分成不同类别,而不同类别需要采用不同的分析方法。因此在进行多因素方差分析之前必须正确判断问题所属类型,否则就可能采用错误的分析方法。
按因素类型进行分类,多因素方差分析可分为固定模型,随机模型及混合模型三类。这几类模型的计算公式基本相同,但其数学模型,假设,统计量,结果的解释等方面均有相当大的差异,我们将在下文中详细介绍,使用时应注意根据实际情况选用适当的模型。
按实验设计分类,多因素方差分析可分为交叉分组和系统分组两大类。这两类计算公式也有些差别,下面我们以两因素方差分析为例,介绍它们试验设计方面的不同点。
交叉分组:实验中,A因素的每个水平都会和B因素的每个水平相遇,因此A,B的地位是完全对称的。这是最常见的实验设计方法。
系统分组:先按A因素的a个水平分为a组,在每一组内再按B的水平细分。一般A因素不同水平的组内B因素的水平可取不同值。例如研究PH值对酶活性的影响,不同的酶可能有不同的最适PH值,因此应对每种酶设置PH值偏高、合适、偏低三个水平,而不同的酶(因素A的不同水平) PH值(因素B)的水平可能是不相同的。
从上面的介绍看出这两种方法适用于不同的问题,必须在实验设计阶段选取适当的方法,才能取得正确的结果。它们的计算方法和公式都是不同的。使用时应加以注意。下面我们具体介绍各种类型的分析方法。
两因素交叉分组方差分析
固定效应模型。首先考虑有重复的情况。线性统计模型为:
xijk=(+(i+(j+((()ij+(ijk, i=1, 2, ……a, j=1, 2, ……b; k=1, 2, ……n
其中:(:总平均值;(i:A因素i水平主效应;(j:B因素j水平主效应;
((()ij:A因素i水平与B因素j水平的交互效应;(ijk:随机误差。
对固定效应模型,应有:
, ,
零假设为:
H01: (i =0, i=1, 2, ……a
H02:βj=0, j=1, 2, ……b
H03:((()ij=0, i=1, 2, ……a, j=1, 2, ……b
备择假设为:
HA: 上述各参数中至少有一个不为0。(这实际上是三个备择假设。)
方差分析的基本思想仍是总变差分解:
即: SST = SSA + SSB + SSAB + SSe
自由度:abn-1 a-1 b-1 (a-1) (b-1) ab(n-1)
均方数学期望分别为:
上述MSA,MSB的均方期望中均不含有交互作用项,这是因为对固定模型来说,交互作用满足:
这说明观测值x只要对i或j中的一个下标求和或求平均,就可以保证交叉项为0。
由于 ,
公式中的x均为平均数,因此上述条件实际保证了在它们的均方期望中不会含有交互作用项。这样,检验两个主效应及一个交互效应的下述三个统计量中,分母全部采用MSe即可。
检验H01,H02,H03的统计量分别为:
(4.17)
, (4.18)
(4.19)
从前述的各均方期望可知,只有当各H0成立时,上述三个分子才是(2的无偏估计量,此时各统计量均服从F分布;若某个H0不成立,则相应的分子将有偏大的趋势,从而使对应的统计量也有偏大的趋势,因此可用F分布上单尾分位数进行检验。
各效应的估计值为:
其中i=1, 2 ……a, j=1, 2, ……b。
实际计算公式为:
(4.20)
(4.21)
(4.22)
(4.23)
或计算:, (4.24)
则:
若使用带统计功能的计算器,可按以下步骤计算:
1°计算排列如下表:
j
i
1 2 ……… b
1
2
a
………
表中最下一行是各列的平均,最右一列是各行的平均。
2°把所有原始数据放在一起,计算样本方差S2,则SST =(abn-1)S2 (4.25)
3°用上表中计算样本方差,则SSST = n(ab-1) (4.26)
4°用上表中计算样本方差,则SSA = bn(a-1) (4.27)
5°用上表中计算样本方差,则SSB = an(b-1) (4.28)
6°SSe = SST - SSST, (4.29)
SSAB = SSST - SSA - SSB (4.30)
完成上述计算后,则可列出以下的方差分析表:
变差来源
平方和
自由度
均方
统计量F
主效应A
主效应B
交互效应AB
误差
总和
把计算所得结果填入上表后,再根据各F统计量的自由度查出其F0.95及F0.99分位数,并将F计算值与相应分位数相比,大于F0.95则在统计量F右上角标一个“*”号;大于F0.99则再加一个“*”号。最后用一句话对上述方差分析的结果加以总结,即哪些主效应或交互效应达到显著或极显著水平,哪些不显著。
如果MSAB小于或约等于MSe,即FAB小于或约等于1,说明此时交互作用不存在,在这种情况下也可把MSAB和MSe合并在一起(即把平方和和自由度都合并)作为σ2的估计量,这样可以提高检验的精确度。具体计算公式如下:
(4.31)
然后可用作统计量FA和FB的分母,对两个主效应进行统计检验(见例题4.7)。注意查表时分母自由度要相应改变。
原料
种类
(A)
温 度(B)
30℃
35℃
40℃
1
41
49
23
25
11
13
25
24
6
22
26
18
2
47
59
50
40
43
38
33
36
8
22
14
18
3
35
53
50
43
38
47
44
55
33
26
29
30
例4.3 为选择最适发酵条件,用三种原料、三种温度进行了实验,得结果如表4.5。请进行统计分析。
表4.5 不同条件下发酵的酒精产量
解:本题中显然温度是一个因素,原料种类是另一个因素。这两个因素各有三个水平。由于它们的影响都是可控制、可重复的,因此都是固定因素。在同样温度、原料下所做的几次实验应视为重复,它们之间的差异是由随机误差所造成的。具体计算过程如下:
用带统计功能的计算器计算。首先计算各处理的平均数,填入下表:
表4.6 各处理平均数表
j
i
1
2
3
1
34.5
18.25
18
23.58
2
49
37.5
15.5
34
3
45.25
46
27
39.42
42.92
33.92
20.12
根据(4.25)~ (4.30)式,有:
把所有原始数据输入计算器,得样本方差S2 = 204.8571,
∴ SST =(36-1)×S2 = 7170.00
把表4.6中间部分9个输入计算器,得样本方差=172.2969
∴ SSST = n(ab-1) = 4×(3×3-1) ×172.2969 = 5513.50
把表4.6中各输入,得样本方差= 64.7575,
∴ SSA = bn(a-1) = 3×4×(3-1)×64.7575 = 1554.18
把表4.6中各输入,得=131.2708
∴ SSB = an(b-1)= 3×4×(3-1)×131.2708 = 3150.50
SSe = SST - SSST = 1656.50
∴ SSAB = SSST - SSA - SSB = 808.82
列成方差分析表,得:
表 4.7 发酵实验方差分析表
变差来源
平方和
自由度
均方
F
原料A
温度B
AB
误差
1554.18
3150.50
808.82
1656.50
2
2
4
27
777.09
1575.25
202.21
61.35
12.67**
25.68**
3.30*
总和
7170.00
35
查F分布表,得:F0.95(2,27)≈F0.95(2,30)=3.316, F0.99(2,27)≈F0.99(2,30)=5.390,
F0.95(4,27)≈F0.95(4,30)=2.690, F0.99(4,27)≈F0.99(4,30)=4.018,
∴FA,FB均达极显著,标上“* *”,FAB只达显著,标上“*”。因此酒精产量不仅与原料和温度的关系极显著,与它们的交互作用也有显著关系。即对不同原料应选用不同的发酵温度。
在固定效应模型中,若各F统计量有达到显著或极显著水平时,常常还需要在各处理间进行多重比较,以选出所需要的条件组合。例如在例4.3中,我们已经发现原料,温度以及它们的交互作用都对酒精的产量有影响,显然我们应进一步找出最优的条件组合以用于生产。这就需要进行多重比较了。如果没有交互作用,可以固定B因素的一个水平,例如取j=1,比较A因素各水平的平均数,得到最优值i*。再固定i,例如仍取为1,比较B因素各水平均值,得到最优值j*。则条件组合A因素i*水平,B因素j*水平就应是所有参加实验的水平组合中最优的。如果有交互作用存在,则一般需要把所有ab个水平组合放在一起比。比较的方法仍与单因素方差分析相同,最常用Duncan法。
例4.4 对例4.3中各处理作多重比较。
解:把各处理平均数从大到小排列(记为x1~x9):
49, 46, 45.25, 37.5, 34.5, 27, 18.25, 18, 15.5
求出各对差值,列成下表:
x9
x8
x7
x6
x5
x4
x3
x2
x1
x2
x3
x4
x5
x6
x7
x8
33.5**
30.5**
29.75**
22**
19**
11.5
2.75
2.5
31**
28**
27.25**
19.5**
16.5**
9
0.25
30.75**
27.75**
27**
19.25**
16.25**
8.75
22**
19**
18.25**
10.5
7.5
14.5*
11.5
10.75
3
11.5
8.5
7.75
3.75
0.75
3
根据公式(4.15),求得:,df=27
查Duncan检验的r值表,求出df=27, k=2~9, α=0.05和α=0.01的r值,并求出临界值R=r, 列成下表:
K
r0.05
R0.05
r0.01
R0.01
2
3
4
5
6
7
8
9
2.91
3.05
3.14
3.21
3.27
3.30
3.34
3.36
11.40
11.94
12.30
12.57
12.81
12.92
13.08
13.16
3.92
4.10
4.20
4.29
4.35
4.40
4.45
4.49
15.35
16.06
16.45
16.80
17.04
17.23
17.43
17.58
将差值表中的数与临界值比较,超过R0.05的标一个“*”号,超过R0.01的标“**”号,一次可核对一条对角线(从左下到右上),因为它们有共同的k值。在第一条最长的对角线上,k=2; 其左上相邻的一条k=3; 余类推,直到左上角最后一个数字,在本题中它的k应取为9。
分析:从这一差值表中可见,x1至x5,除x1至x5外相互间都没有显著差异。但x4,x5与其他3个值差异相对大一些。x6至x9差异均不显著。而x1,x2,x3与x6 ~ x9差异均达极显著。另外,x1,x2,x3以及x7,x8,x9之间的差异都很小。由于现在的数据是发酵产量,显然是越高越好,因此我们主要关心x1,x2,x3。从以上分析中可知,基本上可把x1,x2,x3视为无差异,可在这三组条件组合中,进一步考虑原料成本,原料来源稳定性等其他条件,选一组投入生产。也可对这三组条件增加重复数,进一步检验它们间是否仍有差异。如果实际问题不是要求选最大的数,而是选最小的数,那么根据类似的分析,我们应在x7, x8 , x9对应的三组数中选择。
总之,多重比较的结果分析比较复杂,也比较灵活,需要结合具体数据以及实际问题的要求来进行。这一点请同学们务必注意。
几点注意事项:
1°当交互作用存在时,对固定模型若不设置重复,则无法把SSAB与SSe分开,这样将无法进行任何统计检验。因此在固定模型中有交互作用时,不设置重复的试验是无意义时。
2°对固定模型来说,结论只能适用于参加实验的几个水平,不能任意推广到其他水平上去。
无重复的情况:
刚才我们强调了重复对固定模型方差分析的重要意义,其实重复对所有的方差分析都是相当重要的,这一点我们在后边还会提到。但是重复数每增加1,全部处理的实验就都要多做一次,在工作量方面付出代价也是相当大的。因此,若由经验或专业知识可以断定两因素间确实无交互作用,也可以不设重复,这样可以大大减少工作量。此时线性统计模型变为:
i=1,2,……a, j=1,2, ……b
其中
零假设: H01:(i=0, i=1,2,……a
H02:(j=0, j=1,2, ……b
均方数学期望:
统计量:
其他如结果的解释,计算公式等均与以前一样,只是令n=1即可。
例4.5 在1976-1979四年间四个生产队的小麦亩产量如表4.8所示。各年,各生产队产量是否有显著差异?
表4.8 四个生产队四年小麦田产量(斤)
年 度 (A)
平均()
1976
1977
1978
1979
队
别
(B)
1
546
578
813
815
688
2
600
703
861
854
754.5
3
548
682
815
852
724.25
4
551
690
831
853
731.25
平均
561.25
663.25
830
843.5
724.5
解:本题显然是两因素无重复方差分析。其中生产队和年份各是一个因素。由于生产队对产量的影响主要表现在土地肥沃程度,灌溉条件好坏,耕作习惯差异等方面,在几年内可视为稳定不变的,因此可视为固定因素;而年份对产量的影响则主要体现在气候方面,这是不可重复的,因此应视为随机因素。这样一来,本题实际上成为两因素混合模型方差分析。但由于没有交互效应(这一点最好由专业知识判断,但在本题中专业知识很难判断不同的气候类型对各生产队的影响是否一致,因此我们这里先假设交互作用不存在,后文会提供检验方法),统计计算和检验方法都变得与固定模型完全相同,只是在最后结果的解释上有不同,即固定因素的结果不能推广到其他水平,而随机因素的结果可推广到其他水平。这些差异的原因我们将在随机和混合模型中详细介绍。
先把全部数据输入计算器,得:
∴
再输入各,得
再输入各,得
列成方差分析表:
变差来源
平方和
自由度
均方
F
队别
年度
误差
9111.5
222773.5
5379.0
3
3
9
3037.17
74257.83
597.667
5.082*
124.246**
总和
237264.0
15
查F分布表,得:F0.95(3, 9)=3.863, F0.99(3, 9)=6.992,∴ FA达显著,FB达极显著,分别标以“*”和“**”。即,生产队间产量差异显著,年度间差异极显著。
3. 两因素无重复模型中交互效应的检验。
若由于某种原因不能安排重复,但对是否有交互效应又没有十分把握,则可采用Tukey于1949年提出的一种方法作判断。方法是把残余项(SST-SSA-SSB)再分解,得:
(4.32)
(4.33)
令 若有交互作用,F有偏大的趋势。∴可用上单尾分位数进行检验。
例4.6 判断例4.5中队别与年度间是否有交互作用。
解:
查表,F0.95(1, 8) = 5.32, ∴ 接受H0,可以认为无交互作用。
需要注意的是上述方法虽理论上可行,但在实用中却有很大问题。从(4.32)式可知,SSN的分子实际是两大串数字分别相乘相加再相减,然后再平方。这种计算公式从误差传递的角度看,实在是犯了大忌。因为根据误差传递理论,在相加,相乘过程中,有效数字(即未受误差影响,可以信任的数字)不会增加,而且会集中在头几位非零数字中。而在接下来的相减中,最大的几个非零数字常常是相同的,一减都变成了零,因此有效数字常常会大大减少。在例4.6中,前4位有效数字都损失了,而一般实验中测定的数据有效位数很少有能达到4位以上的。从这一角度说,这种检验方法是非常不可靠的。上述计算只能看作一个计算方法的例子。
综合有关分析,我们可得到以下几点结论:
1°在可能的情况下不采用无重复方差分析;
2°如果必须采用,最好由专业知识保证交互作用不存在;
3°最后没有办法再采用Tukey法进行统计检验,此时应注意计算过程的有效数字位数,尽可能保证结果的可靠性。
4. 无重复方差分析中缺失数据的弥补
方差分析的数据都是按照事先作好的实验设计收集的。但有时由于某种意外的原因,如不可抗拒的自然灾害,实验动物的死亡,操作失误等等,都可能失去一两个实验数据。此时最好的办法当然是重做有关实验来补充,但这有时是办不到的。例如农时一过即不可再种作物,明年气候条件又变化了,无法比较等等。此时如果把整组实验都废弃掉显然是非常可惜的,因此我们需要某种补救的方法。
对于有重复的方差分析来说,丢失一两个数据一般不会造成问题,只要改为按不等重复的方式处理即可。对于无重复的实验设计则必须弥补失去的数据。常用的方法是按照使误差平方和最小的原则来估计缺失的数据。下面以两因素无重复方差分析为例,介绍具体的计算方法。
设缺失的数据为,把它代入SSe的计算公式:
SSe = SST - SSA - SSB
根据最小二乘法,使SSe最小的xij应满足:
0
若用x(i.,x(.j,x(..分别代表去掉未知的xij后的各有关和数,则上式变为:
可解得: (4.34)
上述公式也可从另一思路获得:由线性统计模型有:
xij = ( + (I + (j + (ij,
其中(, (i, (j的估计值分别为:
代入线性统计模型,可得xij的估计值为:
仍用x(i.,x(.j,x(..分别代表去掉未知的xij后的各有关和数,则可得:
这与根据最小二乘法得到的方程是完全一样的,解当然也相同。
若丢失两个数据x, y,仍可采用最小二乘法,令
解上述方程组即可得到x,y的估计值。也可采用迭代法:令,代入(4.34)式,可求出,再把代入(4.34)式,求出y2 ,…… ,这样反复迭代,直到xi-1与xi和yi-1与yi的差很小为止。
几点说明:
1°缺失数据估计出以后,把它填入相应的位置,按一般方差分析的方法计算即可。但自由度会有变化,总自由度应减去缺失的数据个数,SSA;SSB的自由度不变,误差项自由度也相应减小。
2°缺失数据的估计只是一种技术上的处理,它使计算可以进行下去。但是原来的实验数据所应提供的信息却再也找不回来了。因此若缺失数据较多,只好把全部结果报废,勉强分析会得出错误的结论。因此实验时一定要认真,尽量不丢失数据,不能把希望寄托在用计算方法弥补上。
3°弥补的原则是使误差平方和最小,因此处理平方和有偏大的趋势。这相当于引入了一个额外的误差,降低了结论的可靠性。若缺失数据不多,对总的检验结果尚不起太大影响;若缺失数据较多,则应放弃这批数据。
4°在有重复的方差分析中,一般不必进行弥补,只需改用不等重复的计算方法即可。
5.随机效应模型
与固定效应模型相比,线性统计模型本身无变化:
但主效应与交互效应变成了随机变量,它们应满足的条件变为:
因此观察值的方差变为:
。
零假设:
总变差的分解仍同固定模型一样,自由度也不变:
SST = SSA + SSB + SSAB + SSe
df: abn-1 a-1 b-1 (a-1)(b-1) ab(n-1)
均方数学期望变为:
注意上述MSA,MSB的均方期望中,均含有交互作用项,这一点与固定模型是完全不同的。其原因就在于现在是随机模型,交互作用应满足的条件变为。由于现在是随机变量,不再能保证。这样一来,MSA,MSB表达式中均不可能把交互作用项完全消掉,从而也就出现在它们的均方期望中。由于MSA,MSB的均方期望含有交互作用项,检验主效应的统计量也就不能再用MSe做分母,而需要改用MSAB了。
因此,检验各假设的统计量变为:
对检验结果的解释现在也不局限于参加实验的水平,而是可推广到一切水平上。
如果有必要的话,可以根据均方数学期望算出各方差的估计值:
实际计算公式不变,不再重复。对于随机效应模型多重比较是无意义的,因为一般来说处理的效果是无法严格重复的。
与固定模型相同,若FAB的值小于或约等于1,说明交互作用不存在,则可把SSe与SSAB合并。合并方法也与固定模型相同,即为:
(4.31)
然后用作分母构造统计量FA与FB。注意查表时分母自由度也要变为。
6. 混合模型:
不失一般性,我们可假设A因素是固定型,B因素是随机型。线性统计模型仍不变:
xijk = ( + (I + (j + ((()ij + (ijk, (ijk~ NID(0,(2)
条件变为:
但各不是完全独立的,它满足:
即在随机因素的任一水平上均不是独立的。
均方期望:
注意上述均方期望中,固定因素A的均方期望含有交互作用项,而随机因素B反而不含,这跟固定模型和随机模型正好是相反的。造成这种差异的原因还是在满足的条件上:对任意固定,有:,而对固定的i,。这样一来,在MSB的表达式中,和都可保证交互作用被消除掉,从而MSB的均方期望中也就不会有项;但MSA中的却不能使被彻底消去,从而均方期望中也就会出现项。这种均方期望的差异当然会反映在统计量中,即统计量相应变为:
注意上述统计量中由于固定因素的均方期望中有项,要用MSAB作F统计量的分母;而随机因素的均方期望中没有项,要用MSe作F统计量的分母。这正是,而的结果。
固定因素效应估计:
, i=1,2,……a。
方差分量的估计为:
在结果解释方面,固定因素的结论只能适用于参加试验的几个水平,不能推广;而随机因素的结论可推广到它的一切水平上去。其他如变差分解,自由度分解,计算公式,FAB小于或约等于1的处理等均不变,不再重复。
例4.7 为检验三种配合饲料的效果,从三窝仔猪中各选9只,随机分成三组,分别喂以三种饲料。日增重值见表4.9,请对结果作统计分析。
表4.9 仔猪日均增重表
饲料
(A)
窝 别 (B)
1
2
3
1
2
3
1.38
1.30
1.25
1.26
1.23
1.30
1.19
1.23
1.25
1.29
1.32
1.23
1.22
1.28
1.25
1.23
1.18
1.17
1.35
1.40
1.36
1.32
1.28
1.35
1.27
1.31
1.26
解:饲料是固定因素,窝别是随机因素,这是一个两因素交叉分组混合模型。首先把原始数据改写成以下的处理均值
j
i
1
2
3
1
1.31
1.263
1.223
1.266
2
1.28
1.25
1.193
1.241
3
1.37
1.317
1.28
1.322
1.32
1.277
1.232
1.276
1? 把各输入计算器,算得它们的子样方差为根据(4.26)式,
;
2? 把各输入,得其子样方差,根据(4.27)式,得:
;
3? 把各输入,得子样方差,根据(4.28)式,得:
4? 把各原始数据输入,得子样方差S2 = 0.003563,根据(4.25)式,得:
SST = (abn-1)S2 = (3×3×3-1)×0.003563 = 0.09264
5? 由(4.29)式,得:
SSe = SST - SSST = 0.09264 - 0.06636 = 0.02626
6? 由(4.30)式,得:
SSAB = SSST - SSA - SSB = 0.06636 - 0.03116 - 0.03467 = 0.00053
7? 由于a = b = n = 3,各自由度分别为:
dfA = a – 1 = 2
dfB = b – 1 = 2
dfT = abn – 1 = 27 – 1 = 26
dfAB = (a-1)(b-1) = 2×2 = 4
dfe = ab(n-1) = 3×3×2 = 18
8? 把上述计算结果列成方差分析表:
变差来源
平方和
自由度
均方
F
饲料(A)
窝别(B)
AB
误差(e)
0.03116
0.03467
0.00053
0.02626
2
2
4
18
0.01558
0.01734
0.000133
0.00146
117.1**
11.88**
0.091
总和
0.09264
26
查表,得:F0.95(2, 4) = 6.94, F0.99(2, 4) = 18.0
F0.95(2, 18) = 3.55, F0.99(2, 18) = 6.01
F0.95(4, 18) = 2.93
由于FA = 117.1 > F0.99(2, 4), 因此A因素(饲料)主效应达极显著;
由于FB = 11.83 > F0.99(2, 18), 因此B因素(窝别)主效应也达极显著;
由于FAB = 0.091 < F0.95(4, 18), 因此交互效应不显著。
由于FAB <1,为提高检验精度,可将SSAB与SSe合并:由(4.31)式,有:
查表,得:F0.95(2, 22) = 3.44, F0.99(2, 22) = 5.72,
由于FA = 12.77 > F0.99(2, 22), FB = 14.21 > F0.99(2, 22), 因此两因素(饮料与窝别)的主效应均达极显著水平。交互效应显然不显著。
几点注意事项:
1°由于MSAB一般要大于MSe,尤其是交互作用存在时更是显著地偏大,因此若不注意区分是随机因素还是固定因素,就有可能错用统计量,导致错误的结论。因此在两个以上因素的方差分析中,区分因素类型显得更为重要。
2°在随机模型和混和模型中若不设置重复,同样会导致无法把SSAB与SSe分开。此时随机模型仍可对主效应进行检验,混合模型中也可以对固定因素的主效应进行检验。但当交互作用存在时,仅检验主效应是意义不大的,因为很可能是交互作用在起主要作用。因此只要条件容许,不论哪一类模型都应设置重复,除非有可靠的证据证明交互作用不存在。
7. 总结:两因素方差分析表(见表4.10)
表4.10 两因素交叉分组方差分析表
变差
来源
平方和
自由度
固定模型
均方期望
F
A
B
AB
误差
a-1
b-1
(a-1)(b-1)
ab(n-1)
续表4.10
变差
来源
随机模型
混合模型(A固定,B随机)
均方期望
F
均方期望
F
A
B
AB
误差
三、两因素系统分组实验的方差分析。
前面介绍的方法都只适用于交叉分组的实验设计,即A因素的每个水平与B因素的每个水平都会遇到,因此A因素与B因素的地位是完全对称的。但在某些情况下无法采用这样的实验设计。比如进行某种农作物的产量对比实验,A为品种,B为播种期。由于不同品种的最适播期也不一样,采用交叉分组就不太合适,比较理想的方法是根据各自的最适播期分别安排B的水平。这样,先按不同品种分组,然后在每一组内安排自己的播期,这种实验设计方法称为系统分组。其他例如要比较不同菌种的发酵产量,不同酶对同一底物的利用速率等实验中,比较对象对环境条件的要求都是可能有差异的,显然只有让它们各自在自己的最佳条件下工作才能得出正确的结论,因此在这类情况下都需要有系统分组的实验设计方法。
在系统分组实验设计中,首先分组的因素如上述的品种,菌种等称为一级因素,其次分组的(如播期,温度,PH值等)称为二级因素。显然此时两因素不再是对称的,我们的实验目标一般更侧重于测定一级因素的差异。此时的计算方法与分析方法同交叉分组相比均有所不同。为叙述简单,我们下面假定对一级因素A的各个水平,二级因素B的水平数均相同。
线性统计模型:
其中不仅有下标j,还有下标i;表示对于相同的j,不同的i,所代表的二级因素的水平也是不同的。在这里代表二级因素主效应与交互效应之和。由于i不同时二级因素水平j的意义不同,这两个效应已不可能再分开。其他各符号意义同前。
与交叉分组类似,A、B两因素可为固定型,也可为随机型。其应满足的条件与H0也是类似的:
固定型:
,i = 1, 2, …a
,i = 1, 2, …a
, i = 1, 2,…a; j = 1, 2, …b
随机型:
总变差分解为: SST = SSA + SSB + SSe
相应的自由度分解为: abn-1 = (a-1) + a(b-1) + ab(n-1)
这里与交叉分组的不同点是SSB代表B因素的主效应与交互效应之和,已无法再分开。
计算公式为:
(4.35)
(4.36)
(4.37)
(4.38)
将上述各式与交叉分组的(4.20)至(4.24)各式加以比较,即可知SST,SSA的计算公式没有改变,而SSB的(4.37)式其实是交叉分组中的SSST-SSA,因为现在已不需分解B因素的主效应与交互效应。SSe的(4.38)式与交叉分组的(4.24)式相同。
由于系统分组与交叉分组的差别就是前者不需分解B因素的主效应与交互效应,因此采用计算器进行计算时,仍可采用与交叉分组相同的方法计算SST,SSST,SSA,即先计算处理平均数,i水平平均数,然后计算:
1°把所有原始数据放在一起,计算样本方差S2,则SST =(abn-1)S2 (4.25)
2°用处理平均数计算样本方差,则SSST = n(ab-1) (4.26)
3°用i水平平均数计算样本方差,则SSA = bn(a-1) (4.27)
4°令 SSB = SSST - SSA (4.39)
SSe = SST – SSST (4.40)
以下各步骤,如列方差分析表、查表、比较、解释等均与交叉分组相同,不再重复。统计量按以下方法构建:
均方期望及统计量:
对二级因素B来说没有变化:
(4.41)
对一级因素A来说,依B的类型不同而不同:
B固定:
(4.42)
B随机:
(4.43)
上式中,若因素类型为随机型,则和为方差;若因素类型为固定型,则它们都代表平方和,即:
例4.8 比较4种酶在不同温度下的催化效率,特设计如下实验:由于文献记载各酶最适温度分别为30℃,25℃,37℃,40℃,现设定温度水平如下,最适温-5℃,最适温,最适温+5℃。其他条件均保持一致。保温2小时后,测定底物消耗量(毫克)。全部实验重复三次,得结果如下:
温度
酶种类
A1
A2
A3
A4
偏低
适宜
偏高
14.4, 15.2, 13.5
15.9, 15.1, 14.4
13.8, 12.9, 14.6
13.5, 14.4, 15.2
15.1, 16.4, 15.8
15.7, 14.8, 16.0
14.5, 16.3, 15.4
16.4, 18.1, 16.7
15.8, 14.7, 14.1
11.2, 9.8, 10.5
12.5, 10.9, 11.6
10.3, 11.4, 9.9
请进行统计分析
解:由于各种酶的最适温度不同,上述温度水平偏低、适宜、偏高所代表的实际温度是不同的,应采用两因素系统分组方差分析。酶的种类与温度都应为固定因素。酶为一级因素,温度为二级因素。首先计算各平均值,并列成下表:
酶种类
温度偏低
温度适宜
温度偏高
平均()
A1
A2
A3
A4
14.37
14.37
15.40
10.50
15.13
15.77
17.07
11.67
13.77
15.50
14.87
10.53
14.42
15.21
15.78
10.90
首先把各处理平均数,即上表中间的12个数输入计算器,得它们的子样方差为:
由(4.26)式,得:
再把各酶的平均数输入,得子样方差为
= 4.7971
由(4.27)式,得:
再把全部原始数据xijk输入,得子样方差
= 4.6149
由(4.25)式,得:
由(4.39)式:
由(4.40)式:
由于a = 4, b = n = 3, 各自由度分别为:
dfA = a-1 = 3
dfB = a(b-1) = 8
dfe = ab(n-1) = 24
把上述计算结果列成方差分析表:
变差来源
平方和
自由度
均方
F
酶种(A)
129.522
3
43.174
67.64**
温度(B)
16.681
8
2.085
3.266*
误差(e)
15.319
24
0.6383
总和(T)
161.522
35
其中均方 = 平方和 / 自由度,FA = MSA / MSe (4.42式),FB = MSB / MSe (4.41式)。查表,得:
F0.95(3,24) = 3.01,F0.99(3,24) = 4.72
F0.95(8,24) = 2.36,F0.99(8,24) = 3.36
由于FA>> F0.99,因此酶种差异极显著;F0.99> FB > F0.95,因此温度(包括交互作用)造成的差异显著,但未达极显著水平。
如有必要,可对四种酶的四个平均数进行多重比较,也可对同一种酶的三个温度数据进行比较。由于各酶的温度设定不同,对三种温度的总平均数进行比较没有意义。
表4.11 总结:两因素系统分组方差分析表
方差
来源
自由度
固定模型
随机模型
均方期望
F
均方期望
F
A
a-1
B
a(b-1)
误差
ab(n-1)
续表4.11
A固定、B随机
A随机、B固定
均方期望
F
均方期望
F
上表中:
注意:对B因素的检验实际包括主效应和交互效应,它的自由度与交叉分组不同。
四、两个以上因素的方差分析
把两因素方差分析的方法推广到三个或更多个因素理论上不存在问题,但不仅相应的计算过程明显复杂化,更主要的是所需进行的总实验次数也大大增加,因此一般使用较少。当因素多时,实验设计一般改用正交设计的方法,那样可以大大减少实验次数,分析起来也更为方便。正交设计的方法详见实验设计一章。现以三因素交叉分组固定效应模型为例,给出其计算公式及方差分析表。
线性统计模型为:
其中i=1,2……a, j=1,2,……b, k=1,2,……c, l=1,2, ……n.
总变差的分解为:
SST = SSA + SSB + SSC + SSAB + SSAC + SSBC + SSABC + SSe
计算公式和自由度为:
df = (a-1)(b-1)(c-1)
df = abc(n-1)
统计量及均方期望见表4.12。
表4.12三因素交叉分组固定效应方差分析表
变差
来源
平方和
自由度
均方数学期望
F
A
SSA
a-1
B
SSB
b-1
C
SSC
c-1
AB
SSAB
(a-1)(b-1)
BC
SSBC
(b-1)(c-1)
AC
SSAC
(a-1)(c-1)
ABC
SSABC
(a-1)(b-1)(c-1)
误差
SSe
Abc(n-1)
总和
SST
Abcn-1
§4.3 方差分析需要满足的条件
方差分析应满足的条件
要使方差分析达到预期的效果,实验数据必须满足某些先决条件,主要包括以下三点:
1. 可加性。方差分析的每一次观察值都包含了总体平均数、各因素主效应、各因素间的交互效应、随机误差等许多部分,这些组成部分必须以叠加的方式综合起来,即每一个观察值都可视为这些组成部分的累加和。在对每种模型进行讨论前我们都给出了适合这种模型的线性统计模型,这正是可加性的数学表达式。以后的理论分析都是建立在线性统计模型的基础上的,这正说明可加性是方差分析的重要先决条件。在某些情况下,例如数据服从对数正态分布(即数据取对数后才服从正态分布)时,各部分是以连乘的形式综合起来,此时就需要先对原始数据进行对数变换,一方面保证误差服从正态分布,另一方面也可保证数据满足可加性的要求。
2. 正态性。即随机误差ε必须为相互独立的正态随机变量。这也是很重要的条件,如果它不能满足,则均方期望的推导就不能成立,采用F统计量进行检验也就失去了理论基础。如果只是实验材料间有关联,可能影响独立性时,可用随机化的方法破坏其关联性,详见实验设计一章第二节;如果是正态性不能满足,即误差服从其他分布,则应根据误差服从的理论分布采取适当的数据变换,具体方法将在本节后边介绍。
3. 方差齐性。即要求所有处理随机误差的方差都要相等,换句话说不同处理不能影响随机误差的方差。由于随机误差的期望一定为0,这实际是要求随机误差有共同的分布。如果方差齐性条件不能满足也可采用数据变换的方法加以弥补。
条件1的数学表达式是方差分析的线性统计模型,而条件2,3的数学表达式为ε~ NID(0,σ2)。在实用中,条件1,2的满足主要靠理论分析,即如果我们没有理由怀疑数据的正态性,则认为它们是满足的;而条件3则可用一些统计方法进行检验。下面就对具体的检验方法进行介绍。
方差齐性的检验。
在第三章中,我们介绍过两个总体方差是否相等的检验:F检验。但在方差分析中若要对方差齐性进行检验,必然要涉及多个总体方差进行比较的问题。如果一对对进行多次比较,就会像进行多总体均值检验时一样,引起犯第一类错误的可能性大大增高,因此必须采用专门的方法对多个总体的方差一次进行比较。本节中我们介绍三种多总体方差齐性的检验方法,并对它们进行简单比较,同学们可根据需要选用。
对数方差分析。
对数方差分析主要优点是它针对性很强,即只有当各总体方差有差异时才会出现检验通不过的情况;而对其他一些条件,如总体分布是否正态等并不敏感。它的基本思想是把每个要检验的总体即每个不同处理取出的样本再随机地分解成若干子样本,然后分别计算每个子样本的方差并取对数,最后对这些数据进行单因素方差分析。在方差分析中,各处理被视为因素的不同水平,而同一处理的几个子样本的对数方差则被视为重复。由于需要对每个处理的重复观察值都进一步分解成子样本,这种方法要求重复数很多,而这在处理数也较多的情况下是很难实现的,这一点限制了这种方法的应用。对数方差分析的统计假设为:
H0:各处理方差相同;
HA:各处理方差不完全相同。
具体做法为:设共有a个不同处理,每个处理的重复数为ni,则全部观察值可表示为:
xij, i = 1, 2, …a; j = 1, 2, …ni
在对上述的不同处理的样本进一步分割时,应使各子样本的样本含量尽可能接近,且每个处理分割成的子样本组数mi应满足:
(4.44)
各子样本的样本含量记为:
显然应有:
分割后的数据可表示为:
每个子样本的均值和方差分别为:
令 (4.45)
(4.46)
称vij为yij的自由度。然后对yij作方差分析,但要以其自由度vij为权重 。具体公式为:
(4.47)
(4.48)
(4.49)
(4.50)
统计量为:
(4.51)
当H0成立时,上述统计量F服从自由度为的F分布。当H0不成立时,它有偏大的趋势,因此可用分位数对它进行上单尾检验。
总之,对数方差分析方法的优点是比较严谨,针对性也强,检验目标集中在各总体方差是否相等上;缺点是由于要把各样本进一步分为子样本,需要较大的样本容量。
例4.9 用4种方法测定一个沉积的样本中的重金属含量,得结果如下:
方法
测 定 结 果
1
372, 380, 382, 368, 374, 366, 360, 376
2
364, 358, 362, 372, 338, 344, 350, 376, 366, 350
3
348, 351, 362, 372, 344, 352, 360, 362, 366, 354, 342, 358, 348
4
342, 372, 374, 376, 344, 360
这四种测定方法的方差是否相等?
解:各样本平均样本含量为:(8+10+13+6)/ 4 = 9.125
即每个样本大约应分为3个子样本。考虑到各子样本含量应尽量相等,取子样本含量为3。由于原数据应是随机的,分割时不再进行随机化。分组结果如下:
方法
分组结果
mi
1
(372, 380, 382),(368, 374, 366),(360, 376)
3
2
(364, 358, 362),(372, 338, 344),(350, 376, 366, 350)
3
3
(348, 351, 362),(372, 344, 352),(360, 362, 366),(354, 342, 358, 348)
4
4
(342, 372, 374),(376, 344, 360)
2
根据(4.45), (4.46)式计算各组的对数方差(yij)及自由度(vij)
样本
yij
mi
vij
ni
1
3.332, 2.853, 4.852
3
2, 2, 1
8
2
4.862, 4.431, 5.098
3
2, 2, 3
10
3
3.995, 5.338, 2.234, 3.892
4
2, 2, 2, 3
13
4
5.772, 5.545
2
2, 2
6
以各组自由度为权重,求各组平均数:由(4.47)式:
令 (4.52)
则有:v1 = 5, v2 = 7, v3 = 9, v4 = 4
则(4.48)式可改写为:
(4.53)
由(4.53)式,得:
由(4.49)式,得:
由(4.50)式,得:
由(4.51)式,得:
F = 4.906 / 1.624 = 3.021
查表,得F0.95(3,8)= 4.07 > F, 因此接受H0,认为各测量方法的方差相等。
2. 巴勒特(Bartlett)检验
这种方法实际是检验各样本分布的“拖尾”情况是否相同,因此它不仅对各样本方差是否相等敏感,也对各样本是否都服从正态分布敏感。一般来说这是一个缺点,因为当拒绝H0时,我们无法确定是由于方差不全相等引起的,还是由于不全服从正态分布引起的。因此如果我们检验的目标只是各方差是否相等,则应首先检验各总体分布是否均服从正态分布,通过后再做巴勒特检验才比较有把握。但在方差分析中检验方差齐性时,由于我们既需要保证各总体均是正态的(条件2),也需要保证方差齐性(条件3),因此巴勒特检验的这一缺点反而变成了优点。即只要通过了巴勒特检验,正态性和方差齐性就都有了较好的保证,可以不经数据变换直接进行方差分析。反之,若通不过巴勒特检验,则应找出原因并排除,例如排除异常值或进行适当的数据变换。
巴勒特检验的统计假设为:
,(且各总体分布类型相同);
:至少有,(或各总体分布类型不同)。
统计量为: (4.54)
其中为各子样方差以其自由度为权重的加权平均,即:
(4.55)
(4.56)
其他符号意义同前,例如N为总样本含量,α为方差分析的处理数即巴勒特检验的总体数,为各总体样本的子样方差,ni为各总体样本的样本含量。
巴勒特证明了上述统计量K2近似服从χ2分布,其自由度为a-1。从(4.54)式易知,当各相等时,K2=0;当各差异增大时,K2也增大。因此可用χ2分布对K2进行上单尾检验,即当时,拒绝H0。
当各总体样本含量相等时,上述统计量可简化为:
(4.57)
其中 (4.58)
a仍为总体数,即方差分析中的处理数;n为各总体样本共同的样本含量,即方差分析中的重复数,为各子样方差的算术平均数。
注意:当进行巴勒特检验时,一般要求各总体样本含量ni均大于3。
例4.10 调查不同渔场马面鲀体长,结果如下表。请检验方差齐性。
渔场
马面鲀体长(cm)
A
B
C
22.2, 19.1, 20.0, 18.5, 21.4, 19.5
21.6, 22.3, 23.0, 19.2, 20.6, 21.7
17.6, 16.5, 18.7, 19.0, 18.2, 19.4
解:由于各样本含量相等,可使用简化的(4.57),(4.58)式。由所给数据,可算得:
a=3, n=6
由(4.58)式,得:
由(4.57)式,得
K2的自由度为a – 1 = 2,查表,得:,因此接受H0,各渔场马面鲀体长具有方差齐性。
检验
这种方法不如前两种方法严格,它最大的优点是计算简便,只须选取各子样方差中最大的与最小的作一比值,然后再查专门的表格即可。如果只作为方差分析的预备性检验,即检验各处理是否具有方差齐性,它基本上可满足使用要求。
本方法统计量为多个子样方差中最大与最小者的比值,H0为各子样方差相等,HA为至少有一对方差不等。即使在H0成立的条件下,本统计量也不服从任何理论分布,因此必须使用专门编制的临界值表。注意此临界值表与一般F分布表不同,它的相当于普通F分布第一自由度即分子自由度位置的参数是总体数a,而相当于第二自由度即分母自由度位置的则是分子分母自由度中小的一个。具体方法为:
设有取自不同总体的a个子样方差。
令
且记它们的自由度分别为Vmax和Vmin。
则 (4.59)
(4.60)
查专用临界值表(附表16),得。若 > ,则拒绝H0,认为各子样方差不具有方差齐性;否则则接受H0,认为它们具有方差齐性。
例4.11 用法检验例4.9数据的方差齐性。
解:例4.9中,a = 4, n1 = 8, n2 = 10, n3 = 13, n4 = 6
计算可得:
显然:
由(4.59)式,得:
= 233.067 / 54.214 = 4.299
由(4.60)式,得: V = min(5, 7) = 5
查表, = 16.3 >
因此应接受H0,可认为各子样方差相等。
总结:几种检验方差齐性方法的比较。
表4.13 几种检验方差齐性方法的比较
检验方法
优缺点
对数方差分析
巴勒特检验
检验
针对性强,方法严谨,计算较复杂,所需样本量大
除方差齐性外也对偏态敏感,可较好保证正态性及方差齐性。
计算简单,不够严格,需用专门表格。
三、 数据变换
前边曾提到方差分析应满足的三个条件:可加性,正态性,方差齐性。若在这三个条件不满足的情况下进行方差分析,很可能会导致错误的结论。其中第二、第三两条件是互相关联的,因为有些非正态分布,其方差与期望间常有一定的函数关系,如Poisson分布的数据,其期望与方差相等,指数分布的数据,期望的平方等于方差等等。此时显然若均值不等,则方差也不会相等,因此H0不成立时也就不会满足方差分析的条件。在这种情况下,应在进行方差分析之前对数据进行变换,变换主要是针对方差齐性设计的,但对其他两个条件常也可有所改善。由于本课程的特点,我们不介绍变换的数学原理,只介绍常用的变换方法及适用的条件。
1. 平方根变换
用于服从泊松(Poisson)分布的数据。它的方差与均值相等,因此H0不成立时不能满足方差齐性的要求。常见的例子如血球计数值,一定面积内的菌落数,一定体积溶液中的细胞数或细菌数,单位时间内的自发放射数,一定区域内的植物、动物、昆虫数,等等。其特点是每个个体出现在哪里完全是随机的,与其邻居无关。符合这一特点的现象通常服从泊松分布。
方法:把数据换成其平方根,即用代替xij,然后再进行计算。若大多数据值为10左右,个别接近0,可用代替xij。
2. 反正弦变换
用于以百分数形式给出的二项分布数据。即把原二项分布数据乘以100后作为xij,因此数据一般在0~100之间。如果数据集中于30~70之间二项分布本就接近正态分布,此时也可不做变换。但若变化超出上述范围很大则应变换。
方法:令。即先开平方,再取反正弦。也可直接查表得到yij.
变化范围大实际是指p与q相差很大,此时有相当部分观察值大于70或小于30。此时分布是偏的,与正态分布差别很大。若p与q很接近,则数据多在30~70之间,与正态分布差别不大,就可以不变换。
3.对数变换
主要用于指数分布或对数正态分布数据。这些数据的特点是不能取负值,且其标准差σ常与期望μ接近。例如一些描述寿命的数据。
方法:令yij = lg(xij),若大部分数据小于10,个别接近0,可采用yij=lg(xij+1)的变换。然后对yij作方差分析。
4. Box-Cox幂变换。
前三种变换方法都要求我们对总体分布有一种理论上的了解,即知道总体分布的许多特征,从而知道它们服从什么分布。如果对理论分布一无所知,经检验又不是正态分布,则对它的变换常采用幂变换的方法。只要能找到适当的幂值,常常就能成功地将数据正态化。Box-Cox变换就是常用的一种方法。
它的一般形式为:
(4.61)
(4.62)
显然这一方法的关键是确定λ的值。理论证明,使以下对数似然函数L取最大值的λ就是使原始数据正态化的最佳值:
(4.63)
其中n为样本含量,v为自由度。如果xi是一维数据,则v = n - 1;如果是二维数据,则v = n - 2;依此类推。为变换后数据的子样方差,而xi则为原始数据。显然使(4.63)式取最大值的λ不可能用解方程的方法解出,只能用一维搜索计算机程序找出。这是一个典型的优化问题,可使用任何搜索程序对它求解。一般情况下,λ取整数即可。若求出的λ=0,则使用(4.62)式进行变换;若λ不为0,则用(4.61)式进行变换。
需要注意的是并非所有分布形式的数据都可通过数据变换的方法正态化。例如当数据呈双峰状分布(即密度函数有两个峰值)时,就不可能找到一种使它正态化的变换方法。因此变换后的数据仍需要对它是否服从正态分布进行统计检验。
还需要注意一点,即作了变换以后,接着的分析、比较都是对新变量作的。如果希望回到原来的数据上,由于方差、标准差等不能变换回去,因此不能对原数据进行多重比较。
作业
1. 下表是6种溶液及对照的雌激素活度鉴定。指标是小鼠子宫重量。做方差分析,若差异显著则进一步作多重比较。
溶液
鼠
号
对照
I
II
III
IV
V
VI
1
2
3
4
89.9
93.8
88.4
112.6
84.4
116.0
84.0
68.6
64.4
79.8
88.0
69.4
75.2
62.4
62.4
73.8
88.4
90.2
73.2
87.8
56.4
83.2
90.4
85.6
65.6
79.4
65.6
70.2
2.为了调查三块小麦田的出苗情况,在每块麦田中按均匀分布原则设立了一些取样点,每取样点记录30cm垅长的基本苗数。得结果如下表。三块田的出苗情况是否有差异?
田 块
基 本 苗 数
1
2
3
29 24 22 25 30 27 26
20 25 25 23 29 31 24 26 20 21
24 22 28 25 21 26
3. 为选择合理施肥方式,特设计6种施肥方案,各方案施肥成本相同。小区产量如下表。请选择最好的施肥方法。
施肥方案
A
B
C
D
E
F
小区产量
(kg)
12.9, 13.1
12.2, 12.5
13.7, 11.2
14.0, 13.8
15.1, 13.1
14.6, 15.5
12.6, 13.2
13.4, 15.0
12.8, 14.3
10.5, 11.6
10.4, 9.9
12.3, 10.8
14.5, 15.6
17.0, 15.0
16.2, 16.5
15.5, 14.8
13.2, 14.4
13.9, 15.6
4. 随机选取4个小麦品种,施以三种肥料,小区产量如下:
肥料种类
品
种
(NH4)2SO4
NH4NO3
Ca(NO3)2
1
2
3
4
21.1
24.0
14.2
31.5
18.0
22.0
13.3
31.4
19.4
21.7
12.3
27.5
该问题属于哪种模型?从方差分析的结果可得出什么结论?
5. 用两种不同的饲料添加剂A和B,以不同比例搭配饲养大白鼠,每一种饲料添加剂取4个水平,每一处理设两个重复。增重结果如下:(g)
添加剂B
添加
剂A
1
2
3
4
1
2
3
4
32,36
26,24
33,39
39,43
28,22
29,33
30,24
31,35
18,16
27,23
33,37
28,32
23,21
17,19
23,27
36,34
(1)该实验有可能属于哪几种模型?前提是什么?
(2)如果认为是随机模型,设置重复与不设重复对分析结果有无影响?
(3)若实验本身是固定模型,但分析时误认为随机模型,对结论有何影响?若不设重复,又有何影响?
6. 证明无重复两因素固定效应模型中
7. 品种对比实验中希望同时选择适宜的播种量。由于已知的4个参试品种小麦中有一个分孽力强,一个中等,两个偏弱。因此分孽力强的品种选择了播种量为20,25,30斤/亩,分孽力中等的选择了25,30,35斤/亩,分孽力弱的选择30,35,40斤/亩。小区产量如下表。请进行统计分析。
播种量
品种A
品种B
品种C
品种D
20
25
30
35
40
25.4, 20.9, 23.7
28.6, 27.2, 24.8
27.8, 24.2, 23.9
30.4, 29.8, 24.5
32.4, 28.9, 29.1
26.8, 29.4, 29.1
20.2, 24.1, 23.8
21.8, 23.7, 23.5
24.3, 22.6, 20.5
26.4, 30.1, 25.8
29.6, 27.5, 26.4
24.7, 26.1, 30.7
8. 某化工厂用细胞毒理学方法分别检测了接触苯和不接触苯的工人及附近农民的染色体畸变率, 结果如下表。请进行统计分析。
组别
结构畸变(%)
数目畸变(%)
两者均发生(%)
接触工人
不接触工人
农民
5.73
1.40
0.90
9.27
4.60
3.55
2.11
0.19
0.07
9. 某城市从4个排污口取水,经两种不同方法处理后,检测大肠杆菌数量,单位面积内菌落数如下表。请检验它们是否有差别。
排污口
A
B
C
D
处理方法1
处理方法2
9,12,7,5
13,7,10,8
20,14,18,12
17,10,9,15
12,7,6,10
11,5,7,6
23,13,16,21
18,14,19,11