两个总体间的差异如何比较？
研究样本，通过研究样本来分析总体。实际上，所研究的总体往往是无限总体，总体的参数是无法用观察或计算得到的。同理，总体平均数常常无法计算，因而往往用样本平均数作为总体平均数的估计值，因为样本平均数的数学期望等于总体的平均数。

词义解析
离均差是每个观察值的偏离平均数的度量指标。
样本均方是总体方差的无偏估计值。
标准差为方差的正平均根值，用以表示资料的变异度。
抽样分布的标准差又称为标准误，它可以度量抽样分布的变异。

变异系数
标准差和观察值的单位相同，表示一个样本的变异度，若比较两个样本的变异度，则因单位不同或均数不同，不能用标准差进行直接比较。这时可以计算样本的标准差对均数的百分数，称为变异系数。
由于变异系数是由标准差和平均数构成的比数，即受标准差的影响，又受平均数的影响，因此，在使用变异系数表示样本变异程度时，应同时列举平均数和标准差，否则可能引起误解。

正态分布
标准化的正态分布方程就是在正态分布的基础上令U=（y-u)/s,u为正态分布的平均数，s为正态分布的方差。
由于不同的总体的平均数和方差不同，所以将其转换为标准正态分布方程，这样凡要计算一个正态分布的概率只需将y转换为U值，然后查表就可以得出y落入某区间的概率。

假设测验
可从假设的总体里推论其随机抽样平均数的分布，从而可以算出某一样本平均数指定值出现的概率，这样就可以研究样本和总体的关系，从而进行假设测验。 这就是假设测验的基本原理

U测验和t测验
u测验就是根据标准化正态分布的原理进行计算的，u测验是在总体方差为已知，或方差未知单样本容量相当大，可以用样本方差直接作为总体方差进行应用。
同样，t测验也是根据这个原理进行分析的，只不过因为t测验的样本比较小而总体方差又未知，所以就用样本的方差先估算出总体的方差，然后进行分析计算概率的。

成对数据，由于同一配对内两个供试单位的试验条件很是相近，而不同配对间的条件差异又可通过同一配对的差数给予消除，因为可以控制实验误差，具有较高的精确度。

u测验和t测验适用于一个或两个样本平均数的假设测验，方差分析可以用于3个及以上的样本平均数的假设测验。方差分析是用均方来度量试验处理产生的变异和误差引起的变异而已。

方差分析
对一组处理的重复试验数据经对总平方和与总自由度的分解估计出处理间均方和处理内均方（误差均方），并通过F测验处理间所表示出的差异是否真实（比误差大）

方差分析是建立在一定的线性可加模型基础上的，所谓线性可加模型就是指总体每一变量可以按其变异的原因分解成若干个线性组成部分，它是方差分析的理论基础。

方差分析的基本假定

• 可加性
  对于非可加性质料，一般需进行转换，使其效应变为可加性，才能符合方差分析的线性模型
• 正态性
  但是也有研究者发现，数据不服从正态分布对方差分析的结果影响不大
• 误差同质性
  如果发现各处理间的方差相差比较悬殊，一般可用bartlett氏法测验其是否同质；如果不同质，可将方差特别大或变异特殊的处理从全试验中剔除，或者将实验分为几部分，使每一部分具有比较同质的误差方差，以作为较为准确的假设测验。
  卡平方测验可以检验方差同质性

F测验
在一个平均数为u、方差为S的正态总体中，随机抽取两个独立样本，分别求得其均方为s1和s2，将s1和s2的比值定义为F，F值具有s1的自由度和s2的自由度。

在方差分析的体系中，F测验可用于检测某项变异因素的相应或方差是否真实存在，所以在计算F值时，总是将要测验的那一项变异因素的均方作分子，而另一项变异（例如实验误差项）的均方为分母。也就是说如果检测的变异因素存在，那么他的均方就根据自由度的关系而大于限定内的均方。

多重比较
最小显著差数法（实质上就是t测验）、q法、新复极差法LSD

多重比较结果的表示方法
划线法、标记字母法 先将平均数从大到小排列起来，再将不显著的划分为同一组

参数估计法
矩法、最小二乘法、极大似然法

联合方差分析
对用于多年多点实验的分析

相关系数和决定系数
对于坐标点呈直线趋势的两个变数，如果并不需要由X来估计Y,而仅需了解X和Y是否确有相关以及相关的性质（正相关或负相关），则首先应算出表示X和Y相关密切程度及其性质的统计数————相关系数（以r表示相关系数）。决定系数定义为由x不同而引起的y的平方和占y总平方和的比率（用R表述决定系数）
回归系数就是x对y的效应。

偏回归系数
偏回归系数是在其他自变数保持一定时，此变数对因变数的效应。
偏相关系数就是其他变量保持一定是，此变量和因变量的关系。

协变量
通俗的讲，就是在试验过程中对因变量的影响除了自变量外的变量，一些不可控但是能进行测量的变量。在实验设计中，协变量是独立变量，实验者不能操纵，但仍影响实验结果

协方差是在方差分析的基础上，综合回归分析的方法，研究如何调节协变量对因变量的影响效应，从而更加有效地分析实验处理效应的一种统计技术。简单来讲就是对协变量的分析。

回归分析中如果想求得置信区间，可以在进行回归分析时：分析——回归——线性回归——统计——回归系数——误差条形图的表征

协方差分析
直线回归和相关的应用要点（很重要）

偏度
度量数据偏离正态分布的程度，它刻划分布函数对称性，当偏度为正值时，分布向大于平均数方向偏斜，偏度为负值时则向小于平均数方向偏斜；当偏度的绝对值大于2时，分布的偏斜程度严重。
峰度
度量数据服从正态分布时峰的高度，它刻划不同类型的分布的集中和分散程度，当峰度大于3时，分布比较陡峭，峰态明显，即总体变数的分布比较集中。
偏度和峰度是判断正态分布的重要指标

完全随机试验就是简单的单因素方差分析
但是在随机区组试验中，可以用双因素无重复方差分析，因为区组作为局部控制的一项手段，对于减小误差是相当有效的（一般区组间的F测验可以不必进行，因为试验目的不是研究区组效应的）。

条区实验
在多因素实验中由于实施试验的需要，每一因素的各水平都有较大的面积，因而在裂区设计的基础上将同一副处理也连成一片。这样A,B两个因素就互为主副处理，两者的交叉处理为各该水平的处理组合。这就是条区设计。

裂区实验
裂区就是实验因素有主副之分，因此裂区实验的变异的误差项有两个，而一般的随机区组实验误差项只有一个
http://blog.sina.com.cn/s/blog_ab3eddb50102vz3i.html 使用单因素的定制，然后自己设计模型：区组 主效 区组(主效) 副效 主效*副效. 在文件——新建——语法 中进行修改

条区实验
在spss中使用单因素的全因子分析

组内观察值数目相等的单项分组资料的方差分析（spss）：简单的单因素分析
组内观察值数目不相等的单项分组资料的方差分析（spss）：单因素，类型1
组内又分亚组的单项分组资料的方差分析（spss）：单因素，然后将模型修改为{因素 分组(因素) 亚组(因素*分组).}

多因素方差分析中的处理组合间的差异不必管它，

SPSS
许多现实的问题中，仅仅依靠统计描述和简单的统计推断方法是不够的，现实世界中变量间的联系错综复杂，往往要同时考虑多个因素的作用，并为之建立多变量模型。
常用术语
1、因素（Factor）与水平（Level）
因素也被称为因子，就是指可能对因变量有影响的分类变量，而分类变量的不同取值等级（类别）就被称为水平。
2、单元（Cell）
单元也被称为水平组合，或者单元格，是各因素各个水平的组合。
3、元素（Element）
元素指用于测量因变量值的最小单位。根据具体的试验设计，一个单元格内可以有多个元素，也可以只有一个，甚至没有元素。
4、均衡（Balance）
如果在一个试验设计中任意因素个水平在所在单元格中出现的次数相同，且每个单元格内的元素数均相同，则该试验时均衡的；否则，就被称为不均衡。不均衡的试验设计在分析时较为复杂，需要对方差分析模型做特别设置才能得到正确的分析结果。
两个处理的样本量不等，是不平衡试验，不平衡试验用异方差和等方差计算出的t统计量数值是不相同的，而平衡试验用异方差和等方差计算出的t统计量数值是相同的，只是自由度不同，这时两种方法的结果就比较接近，因此实验设计中通常要求做平衡试验。
两个或多个处理下方差相等的情况称为方差齐性，从严格的意义上说，任何两个处理的方差都不会完全相同，我们说方差齐性也只是认为两个处理的方差相差不大，其方差的变异程度不足以影响统计分析结果的正确性，这时采用平衡试验还能够进一步降低方差的差异对统计分析结果的影响。在方差齐性的前提下，平衡试验的统计效率是最高的。如果实验前能够确定方差是非齐性的，则应该对方差大的处理分配较大的样本量。
实际应用中的多数情况方差是齐性的，在实验的处理数目多于两个时，要使用方差分析比较多个处理间平均水平的差异，而方差分析的前提条件是方差齐性，所以等方差的的假设是普遍的。
5、协变量（Covariates）
协变量指对因变量可能影响，需要在分析时对其作用加以控制的连续性变量。实际上，可以简单地把因素和协变量分别理解为分类自变量和连续性自变量。
6、交互作用（Interaction）
如果一个因素的效用大小在另一个因素不同水平下明显不同，则称为两因素间存在交互作用。
7、固定因素（Fixed Factor）与随机因素（Random Factor）
固定因素是指该因素在样本中所有可能的水平都出现了。
随机因素指的是，该因素所有可能的取值在样本中没有都出现，或者不可能都出现。
方差分析模型的适用条件
1、理论上的适用条件
* 各样本的独立性：由于各样本相互独立，来自真正的随机抽样，才能保证变异能够按照模型表达式那样具有可加性（可分解性）；
* 正态性：由于各组的随机误差项被设定为服从正态分布，因此模型要求各单元格的残差必须服从正态分布。
* 方差齐：同样是因为随机误差项，由于在模型中无论何种组合，随机误差项被假定服从相同的正态分布，因此模型要求各单元格都满足方差齐（变异程度相同）的要求。
2、实际操作中对适用条件的把握
（1）单因素方差分析
因模型只有一个因素，设计较为简单，样本有充足的信息量对正态性和方差齐性进行考察，这已经成为标准分析步骤
但是许多人误将正态性理解为因变量应当正态分布，显然这种想法和实际的要求不是一回事。不过由于模型有一定稳健性，只有因变量分布不是明显偏态，分析结果一般都是较稳定的。
至于方差齐性，需要特别指出的是：根据Box的研究结果，在单因素方差分析中，如果各组的例数相同（即均衡），或总体呈正态分布，则方差分析模型对方差略微不齐有一定的耐受性，只要最大与最小方差之比小于3，分析结果是稳定的。
（2）单元格内重复数据的方差分析
以配伍设计方差分析最为典型，此时不需要考虑正态性和方差齐性问题，原因在于正态性和方差齐性的考虑是以单元格为基础单位的，此时每个格子中只有一个元素，当时没法分析了。除配伍设计的方差分析外，交叉设计、正态设计等可以出现无重复数据的情况。但必须指出，这里只有因条件不足，无法考虑适用条件，而不是说可以完全忽视这两个问题，如果根据专业知识认为可能在不同单元格内正态性，方差齐性有问题，则应当避免使用这种无重复数据的设计方案。
当然，从模型的角度讲，实际操作对数据正态性的考虑还有一个办法，就是拟合完毕后作出残差分析图，如果残差呈随机分布，则可知（单元格内）原始数据满足正态条件。
（3）有重复数据的多因素方差分析
由于正态性、方差齐性的考察是以单元格的基本单位，此时单元格数目往往很多，平均每个单元格内的样本粒数实际上比较少。
另一方面，也可能因为只有极个别单元格方差不齐而导致检验不能通过。根据实际经验，实际上在多因素方差分析中，极端值的影响大于方差齐性等问题的影响，因此实际分析中可以直接考察因变量的分布情况，如果数据分布不是明显偏态，不存在极端值，而一般而言方差齐性和正态齐性不会有太大问题，而且也可以基本保证单元格内无极端值。因此在多因素方差分析中，方差齐性往往只限于理论探讨。但对于较重要的研究，则建模后的残差分析时非常重要的。
LSD法：实际上要求将各组均和一个参照水平加以比较。
S-N-K法：两两比较结果则要清楚的多。
1. 首先，它会把各组在表格的纵向上按照均值的大小排序；
2. 其次，在表格的横向各水平被分为了若干个亚组（Subset），不同亚组间的P值小于0.05，而同一亚组各组均数则两两无差异，比较的P值均大于0.05.

当自变量与其他自变量或者协变量相关时，没有明确的方法可以评价自变量对因变量的贡献。例如，含因子A、B和因变量y的双因素不平衡因子设计，有三种效应：A和B的主效应，A和B的交互效应。假设你正使用如下表达式对数据进行建模：
Y ~ A + B + A:B
有三种类型的方法可以分解等式右边各效应对y所解释的方差。
类型Ⅰ（序贯型）
效应根据表达式中先出现的效应做调整。A不做调整，B根据A调整，A:B交互项根据A和B调整。
类型Ⅱ（分层型）
效应根据同水平或低水平的效应做调整。A根据B调整，B依据A调整，A:B交互项同时根据A和B调整。
类型Ⅲ（边界型）
每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:B交互项根据A和B调整。
对平衡实验，那种模型都可以，但是对于非均衡实验，使用类型Ⅰ

R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型Ⅲ方法。
一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，越基础性的变量越应放在表达式前面，因此性别要放在处理方式之前。

方差分析在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的三个水平下都
是来自正态分布的。

QQPlot图是用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布

qq图是正态分位数图，纵坐标是变量的取值，关键是横坐标，参考了以为博友的博客。自己用R写了一个程序验证了一下。基本没问题。
qqplot全名应该是正态分位数图，横坐标的做法：
首先把变量按从小到大的顺序排列，计算变量的长度，即总共有多少个取值，再按顺序计算变量的所有取值的累积百分比，所谓的累积百分比，也就是可以看成是累积概率，比如有10个值，按照从小到大的顺序，第一个值的排序是1， 那么他的所占的百分比就是10%， 紧接着后一个值所占的百分比也会是10%，但是累积概率值为20%， 依次往后计算，因为最后一个值的累积百分比是100%，即等于1，这个值如果计算它的正态分布概率的分位数的话，是无限大的，因此需要对这个值进行修正一下，就是因为这一个值无限大，所以对全体计算出来的累积百分比减去一个适当小的数，修正后的累积百分比与原百分比相差不多，但是回避了最后一个值是1而无法计算的问题。
有了累积百分比之后，相对应的就是累积的概率值。将累积概率值修正后，即得到累积概率，比如以10个值为例，第一个值的累积概率为0.05，查正态分布表，0.05的累积概率，对应的正态分布的Z值为-1.64，这样一次计算，所得的Z值，就是qqplot的横坐标数据。下面以10个数据和30个数据为例说明。

my.qqplot <- function(y){
op <- par(mfrow = c(1, 1))
N <- length(y)
n <- seq(1, N)
xais <- qnorm((n - (.5*N) /N)/ N)

#####中间三句可选，只是为了输出计算过程######
mid <- cbind(sort(y), n, n/N, (n-(.5*N)/N)/N , xais)
colnames(mid) <- c("y", "rank", "cumpercent", "adj-cumper","xaix")
print(mid)
#####中间三句可选，只是为了输出计算过程######

par(mfrow = c(2,1))
qqnorm(y)
plot(sort(y) ~ xais, main = 'my qqplot')
par(op)
}
y <- rnorm(10, mean = 20, s = 50)
my.qqplot(y)
              y rank cumpercentadj-cumper      xaix
[1,]  2.877321   1       0.1      0.05 -1.6448536
[2,]  6.930063   2       0.2      0.15 -1.0364334
[3,] 16.461444   3       0.3      0.25 -0.6744898
[4,] 36.130825   4       0.4      0.35 -0.3853205
[5,] 40.477883   5       0.5      0.45 -0.1256613
[6,] 50.534636   6       0.6      0.55  0.1256613
[7,] 53.425025   7       0.7      0.65  0.3853205
[8,] 54.554269   8       0.8      0.75  0.6744898
[9,]120.496268   9       0.9      0.85  1.0364334
[10,] 125.290253  10       1.0      0.95  1.6448536))

qqnorm(x1)  #数据是否是正态分布的可视化
qqline(x1)
qqplot

(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=nk,labels=1:n,ordered=FALSE)
n是因子的水平个数;k表示每一水平上的重复次数;length=nk表示总观测数;可通过参数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或xA*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值非常小，说明结果非常显著，应该拒绝原假设，认为各因素在不同水平下的试验结果有显著差别，即三种肥料对苹果产量有很大的影响。

zhan<-read.table("zhan.txt",sep="\t",header=T,stringsAsFactors=F)
head(zhan)
install.packages("psych")
library(psych)
zhan_12HS<-zhan[1:243,5:8]
corr.test(zhan_12HS,use="complete")
zhan_13HS<-zhan[244:585,5:8]
corr.test(zhan_13HS,use="complete")
zhan_14HN<-zhan[586:1026,5:8]
corr.test(zhan_14HN,use="complete")
zhan_14SH<-zhan[1027:1311,5:8]
corr.test(zhan_14SH,use="complete")
zhan_15HN<-zhan[1312:1557,5:8]
corr.test(zhan_15HN,use="complete")
zhan_15SC<-zhan[1558:1908,5:8]
corr.test(zhan_15SC,use="complete")