统计学习05：F检验，多组均值比较的假设检验
2020-11-23 本文已影响0人
小贝学生信
简单理解，t检验是两组均值比较的假设检验，P<0.05则表示两组均值存在显著差异。而F检验（一般更常称为方差分析analysis of variance）就是多组均值是否相同的假设检验，P<0.05，则表示多组间均值存在显著差异。
1、方差分析思路
1.1假设检验
（1）零假设：组均值统计量的方差为0，或者说μ1=μ2=μ3=....=μk
（2）计算F统计量：F=组间方差/组内方差；
（3）若F大于临界值，则认为组均值存在显著差异。
1.2变异分解
理解变异分解则可以更好地明白方差分析的思路与缘由。
- 总变异(total variance)：多组数据作为一个整体，计算的离均差平方和；
- 组间方差(between group variance)：组与组之间的方差，一般就用均值代表该组；
- 组内方差(within group variance)：每组内的方差，求和
F statistcs
一般来说：多组数据的总变异都可分为组间方差与组内方差两部分。
而方差分析的目的就是计算：是不是每组间的（均值）差异占总变异的主要部分。考虑两种极端情况
- 组间方差很大，即组均值差别很大；而组内方差很小，即组内数据差别很小，则F值很大（组间方差占总变异的大部分）
2,2,3,3,2;
10,12,11,13,12;
101,100,102,100,99 -
组间方差很小，即组均值差别小；而组内方差大，即每组的数据波动大，则F值小（组内方差占总变异的大部分）
12,25,33,10,30；
15,20,35,31,11
1,60,50,2，34
Variation decomposition
freedom
- 最后补充的一点就是：需要考虑到自由度的影响。组间方差/(组数-1)；组内方差/(总例数-组数)。然后相除，计算F值即可。即F统计量与两个自由度有关，查临界值表时也要注意。
2、实例计算
方差分析有很多种类型，这里演示最常见的两种：单因素方差分析与双因素方差分析
此外R语言里都有现成的函数可直接计算。通过实例计算，对比函数结果，可帮助更好地理解方差分析计算过程。
2.1 单因素方差分析
- 数据来源：R语言中multcomp包提供的Cholesterol reduction for five treatments.即五种治疗方案下胆固醇水平。每组10人，共50个数据。
- 计算所有数据的总变异
sum((cholesterol$response-mean(cholesterol$response))^2)
#[1] 1820.119
- 计算组间方差(组均值与总均值的差值平方，乘以组容量，求和)
aggregate(cholesterol$response, by=list(cholesterol$trt), FUN=mean)
# Group.1 x
#1 1time 5.78197
#2 2times 9.22497
#3 4times 12.37478
#4 drugD 15.36117
#5 drugE 20.94752
mean(cholesterol$response)
#[1] 12.73808
a=aggregate(cholesterol$response, by=list(cholesterol$trt), FUN=mean)
m=mean(cholesterol$response)
sum((a$x-m)^2*10)
#[1] 1351.369
- 计算组内方差(每组数据与本组均值的差值平方，求和)
sum((cholesterol$response-rep(a$x,each=10))^2)
#[1] 468.7504
1351.369 + 468.7504
#[1] 1820.119 即总变异
- 对比函数结果，结果一致
fit0 <- aov(response ~ trt, data=cholesterol)
summary(fit0)
aov()
trt 即为研究的单因素；Residuals 即为组内变异
Df 表示自由度；Sum sq即为变异值；
Mean Sq=Sum Sq/Df ； F value=337.8/10.4；
再计算在自由度(4,45)下，F=337的p值是多少。
2.2 双因素方差分析
即有两个因素（A、B）影响实验结果，分别有不同的水平。根据是否有交互项可分为两类
所谓交互项，即因素A的不同水平可能影响因素B的作用
- 示例数据：来自R的自带数据集 ToothGrowth：对豚鼠进行两种给药（orange juice or ascorbic acid），分别不同的剂量（0.5,1,2）,观察the length of odontoblasts (cells responsible for tooth growth)结果，共60个数据，即特定一种药物，一个剂量重复十次。
- 因素A为不同药物
supp；因素B为不同的给药量dose；分别计算F统计量，得出因素不同水平是否有显著差异。
table(ToothGrowth$supp,ToothGrowth$dose)
# 0.5 1 2
# OJ 10 10 10
# VC 10 10 10
ToothGrowth$dose=factor(ToothGrowth$dose)
（1）无交互项
- 思路：总方差方法不变；组间差异就分为两个（两种药物的两大组，三种药物剂量的三大组）；组内差异就是6组内的方差。
- 总方差
sum((ToothGrowth$len-mean(ToothGrowth$len))^2)
#[1] 3452.209
- 因素A组间方差
aggregate(len, by=list(supp), FUN=mean)
# Group.1 x
#1 OJ 20.66333
#2 VC 16.96333
b=aggregate(len, by=list(supp), FUN=mean)
sum(((b$x-mean(b$x))^2)*30)
#[1] 205.35
- 因素B组间方差
aggregate(len, by=list(dose), FUN=mean)
# Group.1 x
#1 0.5 10.605
#2 1 19.735
#3 2 26.100
aggregate(len, by=list(dose), FUN=mean)
sum(((c$x-mean(c$x))^2)*20)
#[1] 2426.434
- 组内方差：比较复杂。
首先需要当做6组（2*3）考虑。每组分别需要扣除两组因素的影响，再计算与均值的方差
#A因素相较于总均值的影响
supp2mean=rep(b$x[2:1],each=30)-mean(len)
#B因素相较于总均值的影响
dose2mean=rep(c(c$x,c$x),each=10)-mean(len)
#总均值
m=mean(len)
df=data.frame(ToothGrowth$len,supp2mean,dose2mean,m)
sum((df$ToothGrowth.len-df$supp2mean-df$dose2mean-df$m)^2)
#[1] 820.425
- 对比
aov函数结果，结果基本一致
fit1 <- aov(len ~ supp*dose-supp:dose)
summary(fit1)
aov
（2）有交互项
- 如上所述。交互项就是两个因素是有关系的，例如不同剂量的影响受药物差异很大。
- 基本来说，交互项是从上面那个组内方差的进一步分解。（因素A、B的组间方差不变）
- 计算思路：结合上例，计算6组均值，扣除因素A、B单独的效应后再计算方差，是否有显著意义。若有，则可能就是交互效应导致的。
#交互项差异
aggregate(len, by=list(supp,dose), FUN=mean)
# Group.1 Group.2 x
#1 OJ 0.5 13.23
#2 VC 0.5 7.98
#3 OJ 1 22.70
#4 VC 1 16.77
#5 OJ 2 26.06
#6 VC 2 26.14
a=a[order(a$Group.1,decreasing = T),]
df$mean.s.d=rep(a$x,each=10)
sum((df$mean.s.d-df$supp2mean-df$dose2mean-df$m)^2)
#[1] 108.319
#新的组内差异
sum((df$ToothGrowth.len-df$mean.s.d)^2)
#[1] 712.106
- 对比
aov函数结果，结果基本一致
fit <- aov(len ~ supp*dose)
summary(fit)
image.png
- 综上就是结合具体数据进行的方差分析演算。
- 在R中都有直接的函数，但通过自己具体的尝试，可以对返回的结果含义以及F检验理解更深刻。
- 具体的双因素方差分析的数理知识参考ppt，更多种类的R语言方差分析（协方差分析等）R语言操作见之前的笔记
3、两两比较（简单理解）
3.1 原因
- 方差分析的零假设：组均值统计量的方差为0，或者说μ1=μ2=μ3=....
- 所以备择假设是：所有组中，至少有两组不相等（而不是各组都不同）。
- 得到具有显著意义的F统计量之后，有时需要进一步比较是哪些组两两比较存在差异。
3.2 FDR方法
- 最简单的方法就是进行多次t检验，但会导致假阳性率的增加。
例如对一个t检验的假阳性率为5%，对另一个t检验的假阳性率也为5%。那么保证这两次拒绝零假设的概率都正确的概率为：(1-0.05)^2=90.25%<95%，无论p值多小，进行多次的连续推断，总会增加假阳性错误的概率。换一种想法就是对于p=0.05的显著性，对于1000次显著性结果，就会有50个假阳性事件，对于某些分析是不可接受的。因此需要想办法降低50 - 因此就有了常见的FDR值（在生信数据中很常见），可以理解为是p值的校正，一般都会比p值大。结合具体基因差异示例数据
（1）假设某实验对10个基因进行差异分析（10次t检验），得到的10个p值
gene=paste0("gene",1:10)
p.value=round(abs(rnorm(10,mean=0.030,sd=0.1)),4)
df=data.frame(gene,p.value)
raw df
（2）按p值从大到小排序
df=df[order(df$p.value,decreasing = T),]
df$num = 10:1
df
（3）按照Benjamini & Hochberg(BH)法计算FDR值
计算公式为：FDR=p*n/rank。一个注意点就是计算第二个fdr开始，需要与上一个计算值比较，取其中较小的。R代码如下
for (i in 1:10) {
if (i != 1) {
df$fdr[i]=min(df$p.value[i]*(10/df$num[i]),df$fdr[i-1])
} else {
df$fdr[i]= df$p.value[i]*(10/df$num[i])
}
}
df
与R函数
p.adjust结果一致
p.adjust(df$p.value, "BH")
# [1] 0.22810000 0.16522222 0.10925000 0.07528571 0.05766667 0.05766667 0.05766667
# [8] 0.05766667 0.05766667 0.05766667
可以看到，原始p值中6个是小于0.05，但经FDR校正后，全都大于0.05了，即无显著意义。
FDR在生信，尤其海量的基因比较数据中常用。而在医学实验中，进行若干实验组的两两比较，也有很多其它方法。
3.3 其它方法
- 一般来说，Tukey法具有普适性，无论组数多少
- 如果想进行特殊比较，例如对照组与其它实验组的比较，则首选Dunnett法。
这两组方法，在之前R语言笔记中也学到相关函数。
