17实现数据方差分析
方差分析
方差分析是研究一种或者多种因素对实验结果是否有显著性影响的一种有效的统计方法。其中引起观测值不同的原因有很多，不过主要有以下两类：
一是实验过程中随机因素的干扰或观测误差所引起的不可控原因;
二是由于实验中处理方式的不同或实验条件不同引起的可控原因;
主要工作就是将观测数据的总变异按照变异的原因的不同分解为因子效应与实验误差，并对其作出数量分析，比较各种原因在总变异中所占的重要程度，以此作为进一步统计推断的依据。
方差分析中主要的就是 ANOVA 模型的拟合。
有且仅有一个类别型变量，称为单因素方差分析。
当包含两个甚至更多的因子时，便是因素方差分析设计，两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。
注意，表达式中变量的顺序很重要，例如效应会根据表达式中先出现的效应做调整。A 不做调整，B 根据 A 调整，A:B 交互项根据 A 和 B 调整。
假设检验
假设检验是利用样本数据事先在统计假设中按照某种方法进行检验，来判断此假设是否正确。一般地说，对总体某项或某几项作出假设，然后根据样本对假设作出接受或拒绝的判断。使用了一种类似于“反证法”的推理方法，它的特点是： 先假设总体某项假设成立，计算其会导致什么结果产生。若导致不合理现象产生，则拒绝原先的假设。若并不导致不合理的现象产生，则不能拒绝原先假设，从而接受原先假设。但它又不同于一般的反证法。
假设检验的基本思想
假设检验的基本思想是小概率反证法思想。小概率思想是指小概率事件（P < 0.01 或 P < 0.05）在一次试验中基本上不会发生。反证法思想是先提出假设(检验假设 H0)，再用适当的统计方法确定假设成立的可能性大小，如可能性小,则认为假设不成立，若可能性大，则还不能认为假设不成立。
假设检验的检验规则和两类错误
1.确定检验规则 检验过程是比较样本观察结果与总体假设的差异。差异显著，超过了临界点，拒绝 H0；反之，差异不显著，接受 H0。 2.两类错误 （接受或拒绝H0，都可能犯错误） I 类错误——弃真错误，发生的概率为 α （否定了真实的原假设） II 类错误——取伪错误，发生的概率为 β （接受了错误的原假设）
单因素方差分析
银行规定 VIP 客户的月均账户余额达到 100 万元，以此来作为分行的业绩指标。这里的分行就是因子，账户余额是检验的指标，以下数据是从三个分行中随机先后抽取的 7 个 VIP 客户的账户，下面我们通过单因素方差分析来判断三个分行的业绩指标是否有显著性影响。
|分行|账户余额（万元）| |-|-|-|-|-|-| |A1|103\101\ 98\110\105\100\106| |A2|113\107\108\116\114\110\115| |A3| 82 \ 92 \ 84 \ 86 \ 84 \ 90 \ 88 |
> account <- data.frame(X=c(103,101,98,110,105,100,106,113,107,108,116,114,110,115,82,92,84,86,84,90,88), A=factor(rep(1:3, c(7,7,7))))
> account.aov <- aov(account$X~account$A, data=account)
> summary(account.aov)
Df Sum Sq Mean Sq F value Pr(>F)
account$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
DF 表示自由度，sum Sq 表示平方和，mean Sq 表示均方和，F value 表示 F 检验统计量的值（即 F 比），Pr(>F) 表示检验的 p 值，Residuals 为残差。 可以看出，p < 0.05 ，说明有理由拒绝原假设，即认为三个分行之间业绩指标是有显著性差异。
> plot(account$X~account$A)
image.png
在多次重复使用 t 检验的时候会增大犯第一类错误的概率，使得“显著性差异”的结论不一定可靠，所以在进行多次重复比较的时候，需要对 p 值进行调整。
p.adjust.methods
[1] "holm" "hochberg" "hommel" "bonferroni" "BH" "BY" "fdr" "none"
通常在比较次数较多的时候，采用 Bonferroni 方法的效果较好，所以在作多重 t 检验的时候常采用 Bonferroni 方法来进行 p 值的修正。
> pairwise.t.test(account$X, account$A, p.adjust.method="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: account$X and account$A
1 2
2 0.0013 -
3 3.9e-07 6.5e-10
P value adjustment method: bonferroni
p 值在进行修正后会增大，在一定程度上克服了多重 t 检验的缺点。
方差齐性检验是数理统计学中检查不同样本的总体方差是否相同的一种方法。其基本原理是先对总体的特征作出某种假设，然后通过抽样研究的统计推理，对此假设应该被拒绝还是接受作出推断。
常用方法有：Hartley 检验、Bartlett 检验、修正的 Bartlett 检验。
方差分析中有三条前提假设：
不同水平的总体方差相等
因变量是否符合正态分布
待分析的因变量中的个案是否彼此独立
(1) 误差的正态性检验
R 中函数 shapiro.test() 提供 W 统计量和相应的 p 值可用来对数据作正态性检验，当 p 值小于某个显著性水平 α 时，则样本不是来自正态分布的总体，否则来自正态分布的总体。
> attach(account)
> shapiro.test(X[A==1])
Shapiro-Wilk normality test
data: X[A == 1]
W = 0.97777, p-value = 0.948
> shapiro.test(X[A==2])
Shapiro-Wilk normality test
data: X[A == 2]
W = 0.91887, p-value = 0.4607
> shapiro.test(X[A==3])
Shapiro-Wilk normality test
data: X[A == 3]
W = 0.95473, p-value = 0.7724
样本数据在3 个水平下均是正态的。
(2) 不同水平下方差齐性检验（是否相等） 方差齐性检验最常用的方法就是 Bartlett 检验。
> bartlett.test(account$X~account$A, data=account)
Bartlett test of homogeneity of variances
data: account$X by account$A
Bartlett's K-squared = 0.13625, df = 2, p-value = 0.9341
p 值 > 0.05,接受原假设，样本数据是等方差的。
双因素方差分析
考虑了种子的不同品种（A1,A2,A3,A4）以及不同的施肥方式（B1,B2,B3,B4)对种子的产量是否有显著性影响。
||B1|B2|B3| |-|-|-|-|-|-| |A1|325|292|316| |A2|317|310|318| |A3|310|320|318| |A4|330|370|365|
> agriculture <- data.frame(Y=c(325,292,316,317,310,318,310,320,318,330,370,365), A=gl(4,3), B=gl(3,1,12))
>
> agriculture.aov <- aov(Y~A+B, data=agriculture)
> anova(agriculture.aov)
Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
A 3 3824.2 1274.75 5.2262 0.04126 *
B 2 162.5 81.25 0.3331 0.72915
Residuals 6 1463.5 243.92
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根据 p 值可以看出不同品种（因素 A）对产量是有显著性影响的，但是没有理由说明施肥方法（因素B）对产量是有显著性影响的。
以一个检验毒药强弱的实验为例，给 48 只老鼠注射 I、II、III 三种毒药（因素 A），同时有 A 、B 、C 、D 4 种治疗方案（因素 B ），同时每一种因素组合下重复四次测试老鼠的存活时间，分析毒药和治疗方案以及他们的交互作用对老鼠的存活时间是否有显著性影响。(老鼠存活时间（年）)
||A|B|C|D| |-|-|-|-|-|-| |I|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| |II|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| |III|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(time=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), toxicant=gl(3,16,48,labels=c("I","II","III")), cure=gl(4,4,48,labels=c("A","B","C","D")))
>
> rats.aov <- aov(time~toxicant+cure+toxicant:cure, data=rats)
> summary(rats.aov)
Df Sum Sq Mean Sq F value Pr(>F)
toxicant 2 1.0330 0.5165 23.222 3.33e-07 ***
cure 3 0.9212 0.3071 13.806 3.78e-06 ***
toxicant:cure 6 0.2501 0.0417 1.874 0.112
Residuals 36 0.8007 0.0222
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
对于 p 值，因素 toxicant 和 cure 的影响是高显著性的，而交互项对 time的影响不显著。
rats.aov <- aov(time~toxicant+cure+toxicant:cure, data=rats)
summary(rats.aov)
image.png
图中可以再次看出因素 toxicant 和 cure 的影响是显著的。
> shapiro.test(rats$time[rats$toxicant=="I"])
Shapiro-Wilk normality test
data: rats$time[rats$toxicant == "I"]
W = 0.93235, p-value = 0.2656
> shapiro.test(rats$time[rats$toxicant=="II"])
Shapiro-Wilk normality test
data: rats$time[rats$toxicant == "II"]
W = 0.85254, p-value = 0.01485
> shapiro.test(rats$time[rats$toxicant=="III"])
Shapiro-Wilk normality test
data: rats$time[rats$toxicant == "III"]
W = 0.93596, p-value = 0.3024
结果可以看出来因素 toxicant 中“II” 水平不满足正态性要求，其余都满足。
> shapiro.test(rats$time[rats$cure=="A"])
Shapiro-Wilk normality test
data: rats$time[rats$cure == "A"]
W = 0.90078, p-value = 0.1623
> shapiro.test(rats$time[rats$cure=="B"])
Shapiro-Wilk normality test
data: rats$time[rats$cure == "B"]
W = 0.93404, p-value = 0.4249
> shapiro.test(rats$time[rats$cure=="C"])
Shapiro-Wilk normality test
data: rats$time[rats$cure == "C"]
W = 0.88362, p-value = 0.09755
> shapiro.test(rats$time[rats$cure=="D"])
Shapiro-Wilk normality test
data: rats$time[rats$cure == "D"]
W = 0.89749, p-value = 0.1472
结果中看出，因素 cure 均满足正态性要求。
> bartlett.test(time ~ toxicant, data=rats)
Bartlett test of homogeneity of variances
data: time by toxicant
Bartlett's K-squared = 25.88, df = 2, p-value = 2.4e-06
> bartlett.test(time ~ cure, data=rats)
Bartlett test of homogeneity of variances
data: time by cure
Bartlett's K-squared = 13.211, df = 3, p-value = 0.004202
结论看出因素 toxicant 和 cure 均不满足方差齐次性的要求。
