随机微分方程在数理金融中的应用硕士学位论文 精品
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
摘要
复杂数据主要表现在相依、非线性、维数高与不完全观测等,在股市、基因序列和经济等领域中经常出现。
为解决巨型数据集合问题,数据挖掘的理论、方法和技术已应运而生。
而针对诸如怎样同时检验成千上万个基因中哪些基因的表达水平有显著性差异之类的高维统计推断问题,以错误发现率为主要特征的非参数估计方法无疑为其提供了一个有效的解决途径。
本文主要研究考察错误发现率的在各种参数模型和非参数模型下的控制检验方法,全文共分为四章。
文章首先介绍了所选取课题的背景和意义,以及国内外在该方向的研究现状。
在多重假设检验的背景下,给出了错误发现率的定义,提出利用p值进行假设检验,并在假设检验独立和相依的情形下对错误发现率的控制方法进行了探讨。
在研究错误发现率的控制方法时,发现在处理多重假设检验问题时,核心的问题是如何估计真实零假设的个数,因此本文采用经验贝叶斯估计来估计它的值。
在参数混合模型和非参数混合模型中研究真实零假设的估计问题是本文的核心内容。
针对正态混合分布模型和Beta混合分布模型两种参数混合模型,文章采用矩估计方法和基于p值的最小二乘估计方法进行研究;在研究非参数混合模型时,分别介绍了最小二乘估计方法、Beta分布拟合模型和Beinstein 多项式拟合模型的方法。
文章的最后以Hedenfalk报告的一组乳腺癌患者的基因数据为例进行仿真研究,发现错误发现率为微阵列数据的多重假设检验提供了合适的错误控制指标。
关键词:错误发现率;多重假设检验;p值;非参数估计;微阵列数据
- I -
Abstract
Complex data always appear in the stock market, gene sequences, economic and other fields, which mainly show the characteristic of dependent, nonlinear, high dimension and incomplete observations. In order to solve the problem of huge data collection, the theories, methods and techniques of data mining are proposed. While how to examine the high-dimensional statistical inference problem, such as the significant differences of expression levels in thousands of genes, the non-parametric estimation of false discovery rate provide an effective solution.
This paper mainly investigate the test method based on the false discovery rate of various parametric model and non-parametric model, which is divided into four chapters. Firstly, this paper introduce the background and significance of the topic, and the current studies in this direction at home and abroad. Under the background of multiple hypotheses testing, the paper describe the definition of the false discovery rate, propose using the p-value to test the hypothesis testing, and discuss the controlling method of the false discovery rate when the hypotheses testing is independent or dependent. When we investigate the controlling method of the false discovery rate and studied the multiple hypothesis testing problem, we find that the central problem is how to estimate the number of true null hypothesis, so this paper use the empirical Bayes estimation to estimate its value. Investigating the estimation of true null hypothesis in the mixing parametric model and non-parametric model is core of the dissertation. Aiming at the mixed normal distribution model and Beta mixture distribution model, This paper use the method of moment estimation and least squares estimation method based on the p-value to estimate its value; On studying the
non-parametric mixture model, the paper introduce the least square estimation method, Beta distribution fitting model method and the Beinstein polynomial fitting model method. Finally, the paper conduct the simulation research based on a group of patients with breast cancer gene data by Hedenfalk, and find that the false discovery rate is able to provide a suitable error control targets for the multiple hypothesis testing of microarray data.
Keywords: false discovery rate, multiple hypotheses testing, p-value, non-parametric estimation, microarray data
- II -
- III -
目 录
摘 要 ..................................................................................................................... I Abstract ................................................................................................................... I I
第1章 绪 论 (1)
1.1 课题研究的背景及意义 (1)
1.2 国内外在该方向的研究现状 (1)
1.2.1 国外对错误发现率的研究现状 (1)
1.2.2 国内研究现状 (3)
1.3 本文拟研究的主要内容 (3)
1.4 创新点 (3)
第2章 错误发现率的多重检验方法 (5)
2.1 多重假设检验的错误测度 (5)
2.2 P 值的定义、性质和计算方法 (6)
2.3 独立情形下基于FDR 控制的检验方法 (7)
2.4 相依情形下基于FDR 控制的检验方法 (8)
2.5 真实零假设的个数0m 或比值0π的估计 (9)
2.5.1 -λ估计 (9)
2.5.2 经验贝叶斯估计 (11)
2.6 本章小结 (12)
第3章 参数混合模型和非参数混合模型的估计 (13)
3.1 引言 (13)
3.2 正态分布混合模型 (13)
3.3 Beta 分布混合模型 (17)
3.4 非参数混合模型的估计 (21)
3.4.1 最小二乘估计 (22)
3.4.2 Beta 分布拟合模型 (23)
3.4.3 Beinstein 多项式拟合模型 (25)
3.5 本章小结 (26)
第4章 错误发现率的估计方法的应用 (27)
4.1 引言 (27)
4.2 微阵列数据实例研究 (27)
4.3 本章小结 (28)
结论 (30)
参考文献 (31)
哈尔滨工业大学学位论文原创性声明和使用权限 .......... 错误!未定义书签。
致谢 .................................................................................. 错误!未定义书签。
- IV -
第1章绪论
1.1 课题研究的背景及意义
复杂数据主要表现在相依、维数高、非线性与不完全观测等,经常出现在股市、基因序列和经济等领域中。
在研究处理低维的简单数据时,采用传统的数理统计方法是有效的,但在研究比较复杂的数据时,就会变得比较困难。
因此,“复杂数据的统计推断问题”已被列为我国统计学研究的重点课题。
随着科学技术的不断发展,在实际的统计研究过程中,出现了越来越多的大型数据集合问题。
在研究巨型数据的高维统计推断问题时,以错误发现率为主要特征的非参数估计方法为其提供了一个有效地解决途径。
在巨型数据问题的统计分析中,错误发现率( false discovery rate, FDR)有着非常重要的作用,现已被越来越多地应用在微阵列(Microarray)数据研究和功能磁共振成像(Functional magnetic resonance imaging, fMRI)等领域。
以微阵列数据研究和功能磁共振成像(fMRI)为代表的现代生物技术已经给医学界的研究带来了很大的影响。
由于错误发现率可以为大规模数据多重检验中的错误控制提供一个合适的测量标准,因此在微阵列数据的研究中,研究者通常采用错误发现率(FDR)来控制多重假设检验的错误率。
例如在研究基因表达的差异性试验中,假设我们挑选了R个差异表达的基因,其中有S个是真正有差异表达的,另外有V 个其实是没有差异表达的,也就是说是假阳性的。
在试验中我们希望错误比例
α),在统计学意义上,这就等价
=
.0
V/不能超过某个预先设定的值(比如05
R
于控制FDR不能超过%
5.
1.2 国内外在该方向的研究现状
1.2.1 国外对错误发现率的研究现状
多重假设检验的统计显著性问题已经引起了许多统计学者的注意。
1995年,Benjamini和Hochberg在研究多重假设检验时首次提出了错误发现率的概念,并在多重检验中对它的控制方法做了研究,给出了计算方法[1]。
然而,由于当时没有学者研究大规模数据,因此并未受到重视,甚至还受到广大学者的质疑。
若干年后,随着微阵列数据研究的不断发展,大规模数据的频繁出现使得FDR有了实际的应用,错误发现率的理论和应用研究也在逐渐走向成熟。
FDR(false discovery rate)的定义如下:
⎪⎩
⎪⎨⎧=≠⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛+=0,00,R R R V E S V V E FDR 上式中的V 和S 分别表示m 个假设检验中错误拒绝和正确拒绝检验的个数,R 表示m 个假设检验中总的拒绝原假设的个数,()⋅E 表示数学期望。
Storey 和Tibshirani(2003)提出了阳性错误发现率( positive false discovery rate ,pFDR)的定义,并在DNA 微阵列数据试验应用过程中,分别给出了统计数据独立和相关条件下的程序计算过程[2]。
pFDR 的定义为:
⎥⎦
⎤⎢⎣⎡>=0|R R V E pFDR 其中V 和R 的含义与上文相同。
比较FDR 和pFDR 两者的定义可知,pFDR 是FDR 的一种特例。
设假设检验H 的检验统计量为T ,分别假设0=H 和1=H ,令()00Pr π==H 和()()1,0,|Pr ===≤j t F j H t T j 分别表示检验统计量T 的零分布和相间分布。
同时进行m 次试验。
也就是说,考察m 个假设检验:m H H ,⋯,1及其检验统计量m t t ,,1⋯. 对每个i ,分别假设0=i H 和1=i H . 假定对每个i ,都有()00Pr π==i H 和()()1,0,|Pr ===≤j t F j H t T j . m t t ,,1⋯被当做T 的一个样本,且具有混合分布 ()()()t F t F t F 1000)1(ππ-+= (1-1)
设全体试验的拒绝域的集合为Γ。
未被发现的错误率( false non-discovery rate ,FNR)首次被Genovese 和Wasserman(2002)[3]提出。
从参考文献[4]和[5]中,我们可以得到正错误发现率(pFDR)和未被发现的错误率(FNR)的贝叶斯解释:
()()
()
Γ∈Γ∈=Γ∈==T T T H P pFDR F F Pr Pr |000π ()()
()Γ∉Γ∉-=Γ∉==T T T H P pFNR F F Pr Pr 1|100π 上式中的分母()Γ∈T F 0Pr 和()Γ∉T F Pr 可以由经验分布估计得出结果,有时也会从已知的或者由采样的方法得到的零分布中得到结果。
如果0π可以由检验统计量m t t ,,1⋯估计,那么pFDR 和pFNR 就是可以估计的。
Allison 等人(2002)采用有限Beta 混合模型,利用这些数量模拟了微阵列数据分析中的P 值[6]。
关于多重假设检验问题的研究,也受到了国际著名统计学家的高度重视,且已编入了国际统计学的教材中。
Erich Lehmann 编著的《Theory of Point Estimation 》
和《Testing Statistical Hypotheses 》是世界各国培养统计学研究生的标准教材,被世界各国的大学广泛采用。
2005年,Lehmann 还撰文提出了k-族错误率(k-FWER)的概念。
另外,斯坦福大学统计系教授Bradley Efron 也对此问题作出了深入的研究,并在许多重要报告中介绍了FDR 的应用成果[7-9]。
1.2.2 国内研究现状
在国内统计学研究中,目前对多重假设检验中错误发现率问题的研究才刚刚起步。
黄丽萍等(2003)以脑功能磁共振成像(fMRI)为实验,对多重假设检验的FDR 控制方法进行了研究,他们利用计算机编程技术对FDR 控制方法进行了详细的研究,并在功能磁共振成像(fMRI)数据分析中加以应用[10]。
缪柏其(2005)和朱钰(2005)介绍了FDR 控制检验方法取得的显著成果[11]。
东北师范大学郭建华教授指导的裴艳波(2005)的硕士论文对多重假设检验问题中关于三种错误测度--FWER,FDR 和pFDR 及其控制方法进行了较全面的介绍[12]。
此外,苟鹏程(2006)对微阵列数据的多重比较进行了探讨[13]。
1.3 本文拟研究的主要内容
本文主要研究错误发现率的非参数估计方法,并以微阵列数据为实例进行仿真研究。
在第二章中,我们从多重假设检验的错误测度的角度出发,引入错误发现率的概率意义,研究了p 值的定义和性质,并着重介绍真实零假设的个数0m 或比值0π的估计方法;在第三章,我们详细介绍比值0π在参数混合模型与非参数混合模型下的估计方法;第四章以微阵列数据为例,进行仿真研究,并得出相关结论。
1.4 创新点
本文的创新点在于:
首先,本文在多重假设检验的背景下,介绍了错误发现率的定义,并提出利用p 值进行假设检验;
其次,在对正态混合分布模型和Beta 混合分布模型两种参数混合模型进行研究时,文章采用矩估计方法和基于p 值的最小二乘估计方法进行研究;在研究非参数混合模型时,分别采用最小二乘估计方法、Beta 分布拟合模型和Beinstein 多项式拟合模型的方法进行研究;
最后,在以Hedenfalk 的乳腺癌微阵列数据作为实例对微阵列数据进行仿真研究时,本文采用置换检验的方法对错误发现率的控制方法进行研究,得到合理
的实验结果。
第2章 错误发现率的多重检验方法
2.1多重假设检验的错误测度
在研究多重假设检验问题时,最核心的内容就是如何控制总体检验所犯的错误。
由于涉及多重检验,因此情况将变得非常复杂。
例如,同时对m 个假设进行检验,分别记为m i H i ,⋯=,1,. 如果原假设0i H 为真,则令0=i H ,否则令1=i H . 记{}0:0,i i H M =={}1:1i i H M ==. )(#),(#1100.M m M m ==, 即01,m m 分别为0M 和1M 中含有的元素的个数。
显然有10m m m +=. 对于这m 个检验结果的分类见表2-1.
其中,R 表示拒绝总数,即m 个检验中显著性假设的个数,是一个可观测的随机变量;V 表示m 个检验中犯第Ⅰ类错误的个数;T 表示犯第Ⅱ类错误(假阴性)的总数,T 和V 均为不可观测的随机变量。
在实际的检验过程中我们发现,表2-1中的一些量,例如0,,,,m V S U T 是不
可观测的。
在多重假设检验中,为了衡量检验总体的第Ⅰ类错误,我们必须首先要寻找一种比较合理的错误测度,然后进一步研究该错误测度的控制检验方法,以达到尽可能多地发现显著性假设的目的。
这里我们主要介绍错误发现率(FDR) 的定义及其衍生出来的各种相关测度。
定义2.1:FDR }{0,V E I R R ⎛⎫≡> ⎪⎝⎭
称为错误发现率(False discovery rate)。
1995年,Cahgeton 和Peshereg 提出了错误发现率的概念。
下面是由错误发现率衍生出来的各种相关概念。
定义2.2:cFDR(r),V E R r R ⎛⎫≡= ⎪⎝⎭
称为条件错误发现率(conditional FDR). 定义2.3:eFDR(r),EV r ≡
称之为经验FDR(empirical FDR). 定义2.4:mFDR ,EV ER
≡称之为边缘FDR(marginal FDR). 定义2.5:pFDR 0,V E R R ⎛⎫≡> ⎪⎝⎭
称之为阳性FDR(positive FDR). 定义2.6:FNR }{0,T E I W W ⎛⎫≡> ⎪⎝⎭
称之为假非发现错误率(False non-discovery rate).
定义2.7:pFNR 0,T E W W ⎛⎫≡> ⎪⎝⎭
称之为阳性FNR(positive FNR). 这些衍生的错误测度与FDR 之间的关系可有下列式子表示出:
()Pr 0,FDR pFDR R =⋅>
()()()11
Pr Pr m
m r r V FDR E R r R r cFDR r R r R ==⎛⎫==⋅==⋅= ⎪⎝⎭∑∑ ()()()11Pr 0Pr 0m m r r V pFDR E R r R r R cFDR r R r R R ==⎛⎫==⋅=>=⋅=> ⎪⎝⎭∑∑ 且当0m m =时,有1cFDR mFDR pFDR ===和
()()Pr 0Pr 0 1.FDR R V FWER =>=>=≤
2.2 P 值的定义、性质和计算方法
为了能够直观的得到接受或拒绝原假设的置信程度,我们通常采用P 值来研究。
在多重假设检验的研究中,采用P 值进行假设检验已经成为国际上比较流行的方法。
因此,在研究模型的估计方法之前,有必要先研究下P 值的定义和性质。
定义2.8:设检验统计量为X ,样本观测值为x ,对于一族拒绝域}{,Γ统计量X x =的P 值可以定义为:
()}{()}
{0:min Pr .x p x X H Γ∈Γ≡∈Γ 在实际的假设检验中,由定义2.8所得到的P 值,如果0.05p <, 说明检验结果是显著的;如果0.01p <, 则说明检验结果非常显著。
下面给出P 值的计算方法和作用,并不加证明的给出P 值的性质。
(1)P 值的计算方法
当0H 为真时,统计量X 的值x 可由样本数据计算出,根据检验统计量X 的实际分布,可以求出P 值()x p . 具体地讲,就是:
1.左侧检验的P 值是统计量X 小于样本统计值x 的概率,即:
()}{
0Pr p x X x H =<; 2.右侧检验的P 值是统计量X 大于样本统计值x 的概率,即:
()}{0Pr p x X x H =>;
3.当统计量的分布具有对称的性质时(例如正态分布,t 分布等),双侧检验的P 值是统计量X 落在样本统计值x 为端点的尾部区域内的概率的2倍,也就是说:当x 位于分布曲线的右侧时,有
()}{}{
00Pr 2Pr ,p x X x H X x H =>=>
当x 位于分布曲线的左侧时,有
()}{}{
00Pr 2Pr .p x X x H X x H =>=<
(2)P 值的性质
1.如果原假设为真,那么由定义
2.8计算出的P 值满足()0,1区间上的平均分布,即()0,1P U ;
2.如果原假设非真,则P 值的分布不易确定,但由P 值的统计意义可知,其分布具有递减的趋势。
(3)P 值的作用
在假设检验中,我们先利用样本数据计算出P 值,然后将P 值()x p 与提前给出的检验水平α比较,得出检验的结论:
1.如果(),p x α≤则在显著水平α下接受原假设;
2.如果(),p x α>则在显著水平α下拒绝原假设。
在实际实验过程中,若()p x α=,则可以适当提高样本的大小,再次进行检验。
2.3 独立情形下基于FDR 控制的检验方法
在多重假设检验中,利用P 值来研究错误发现率的控制方法有很多。
在这一节,我们先研究独立情形下基于FDR 控制的检验方法。
Benjamini 和Hochberg 在提出错误发现率的概念的同时,给出了FDR 最初的检验方法,记为BH 法。
BH 法:
设m 个假设检验1,m H H …,对应的P 值分别为1,,m p p …,将他们从小到大排序,得到()()1,,m p p …,其中()i p 对应于.i H 对于给定的检验水平q ,令
()max :,i i k i p q m ⎧⎫
=≤⎨⎬⎭⎩
则拒绝1,k H H …,对应的原假设。
实际上,当时提出的BH 法只是用来控制总体的错误测度(FWER)。
由下面的定理我们可以发现,如果检验水平q 已知,那么该方法就可以有效地控制FDR 。
定理2.1 :【Benjamini and Hochberg (1995)】[14]如果多重假设检验的统计量所对应的P 值相互独立,且具有连续的分布,q 为给定的检验水平,那么BH 法
控制q q m
m FDR ≤≤
. 受BH 法的启示,Benjamini 和Liu(1999)提出了一个step-up 的错误发现率的检验方法,记为BL1法[15]。
BL1法:
取1
111min 1,,1,1m i i m q i m m i δ-+⎡⎤⎛⎫⎢⎥=--≤≤ ⎪⎢⎥-+⎝⎭⎣⎦
经过计算可以知道,i δ是单调上升的,即10 1.m δδ<≤≤≤…令()}
{
min :,i i k i p δ=> 则拒绝()()11,k H H -…,所对应的零假设。
注:在上面的BL1方法中,如果不存在这样的k , 那么拒绝所有的原假设。
定理2.2 :【Benjamini and Liu(1999a)】如果多重假设检验的统计量所对应的P 值相互独立,且具有连续的分布,,则BL1法控制,FDR q ≤ 其中q 为提前给定的检验水平。
由定理2.2可知,在相互独立的条件下,BH 法把FDR 的水平控制在0/.m q m 若0m 已知,则可令'0/q q m m =⋅取代BH 法中的检验水平q ,从而可以更精确地控制FDR 在q 水平内。
2.4相依情形下基于FDR 控制的检验方法
在上一节,我们讨论了独立情形下FDR 控制的检验方法,但是在实际的试验过程中,统计量一般会具有着这样或那样的依存关系,从而使得上面研究的控制方法是无效的。
因此本节将介绍在统计量对应的P 值相依的条件下FDR 的控制方法。
针对多重检验中检验统计量自由分布的情形,我们有下述检验方法,由于该方法是由Benjamini and Liu 提出来的,因此叫做BL2法[16]。
BL2法:
令(){
}2
min 1,/1i d m q m i =⋅-+, 有101m d <≤≤≤…d , 令(){}
min :i i k i p d =>,那么拒绝()()11,k H H -…, 对应的原假设;若不存在上述条件的k , 则拒绝所有原假设。
定理2.3 :【Benjamini and Liu(1999b)】上述针对分布自由的检验统计量的BL2法,有q FDR ≤.
针对多重检验中检验统计量自由分布的情形,还有下述的FDR 控制方法,
该方法由Benjamini and Yekutieli 提出,因此记为BY 法。
BY 法:
令()1min :1i m j i k i p q m j =⎧⎫⎪⎪
⎪⎪
=≤⎨⎬⎛⎫⎪⎪ ⎪
⎪⎪
⎝⎭⎩⎭∑ ,则拒绝()()11,k H H -…, 所对应的原假设。
注:在上述检验方法中,如果不存在这样的k ,则不拒绝任何原假设。
定理 2.4:【Benjamini and Yekutieli(2001)】上述对于多重检验自由分布的step-down 的FDR 控制方法控制FDR 在0/m q m 水平[17]。
2.5 真实零假设的个数0m 或比值0π的估计
通过上文在独立情形和相依情形下基于FDR 控制的检验方法的研究,我们可以知道,在多重假设检验中,如果真实零假设的个数0m 或者比值m m /00=π已知,那么就可以根据检验统计量之间相依或者独立的关系,采用上文介绍的检验方法来控制FDR. 然而在实际研究中,0m 或者0π往往是未知的,因此,最重要的问题就是如何估计0m 的值,或者等价的估计0π的值。
本节就来研究这个问题,我们分两种方法进行具体研究。
2.5.1 -λ估计
基于P 值在不同假设条件下的分布差异性,Storey(2002)提出了一种0π的估计方法,记为-λ估计方法[18]。
若假设m H H ,,1⋯同分布,设m p p ,,1⋯为m 个假设m H H ,,1⋯所对应的P 值。
对()1,0∈∀λ,我们记()}{λλ>=i p W #,那么0π可由下式估计出:
()()
.10λλπ-=
∧
m W (2-1)
由上式可以看出,λ的取值不同,由(2-1)式所得到的0π的估计值就不同,且所得到的估计值都比真实值偏大,这是因为()1,0∈∀λ,有
()[]()∑=>=m
i i p W E 1Pr λλ
()()[]∑==>+=>=m
i i i i i H p H p 1
1,Pr 0,Pr λλ
()()()[]∑==>-+=>=m
i i i i i H p H p 1
001Pr 10Pr λπλπ
()()()1Pr 1100=>-+-=H p m m λπλπ, 从而有
()()1Pr 110
00=>--+
=∧
H p E λλππλπ ()0
1
1
0011π
λ
ππλ≥--+
=⎰dp p g ,
上式中的()p g 1表示备择假设下P 值的密度函数。
由P 值的性质可知,密度函数()p g 1是渐进递减的,而且显然有()01≥p g . 因此由上式可以看出,当λ减小
时,0π的误差
()dp p g ⎰--1
1
011λλ
π就会变小,反之则变大。
而且由()()
λλπ
-=∧1/0
m W 可以看出,当λ增大时,∧
0π的方差就会增大,这就造成了估计值∧
0π的不稳定性。
那么如何才能寻找一个合适的λ,使得估计值∧
0π达到最优呢?2002年,统计学家Storey 提出了一个选取λ的计算方法:考虑使均方误差
()()200⎩⎨⎧⎭⎬⎫-=∧πλπλE M S E (2-2)
最小化的λ取值。
由于上式中的0π未知,我们可以用()λππλ
∧
∧=00min p
取代(2-2)
式的0π,这是因为对()1,0∈∀λ,估计值()λπ∧
0都偏大,于是有
()()2
10*01∑=∧∧∧⎪⎩
⎪⎨⎧⎪⎭⎪
⎬⎫-=K
k p k K MSE πλπλ
其中()λπk
*
∧表示第k 次对P 值样本进行抽样后,采用(2-1)式重新计算得到的估计
值。
从而最优λ为
(),min arg ⎩
⎨⎧⎭⎬⎫=∧
∧
λλλMSE 从而可以得到最优-λ估计
.00
⎪⎭
⎫
⎝⎛=∧∧
∧λππb
由于()1,0∈λ,因此我们可以考虑采用格点法,即在()1,0区间上等距离地抽取有限个λ值,然后利用(2-1)式计算最小化均方误差[19]。
2.5.2 经验贝叶斯估计
在对微阵列数据进行研究时,Efron, B. and Tibshirani, R. (2002)[20]提出可采用经验贝叶斯方法来估计FDR. 令()0Pr 0=≡H π表示不同条件下基因表达无差别的概率,则()0111Pr ππ-==≡H 表示基因表达存在差别的概率。
我们采用10,f f 来表示零假设和备择假设检验下检验统计量X 的密度函数,对应的分布函数分别为10,F F . 则检验统计量X 的密度函数可以表示为
()()().1100x f x f x f ππ+=
计算后验概率,有
()()()
()
,1Pr 111x f x f x X H x p π=
==≡
()()()
()
.0Pr 000x f x f x X H x p π=
==≡
如果()x f 是已知的,或者已经被估计出来,记为()x f ∧
,则由
()()
()
R x x f x f x p ∈∀≥-
=∧
∧
,01001π
得到不等式
()()
.min
00x f x f x
∧
≤π 从而得到0π的一个估计式
()()
.min
00x f x f x
∧
∧
=π 上式也可以改写为
()()
,min
0*
x F x F m x
∧
∧
∧=π
其中()x F m ∧和()x F ∧
0为对应的经验分布函数。
2.6 本章小结
在第一节中,我们介绍了多重假设检验中错误测度的定义,给出了错误发现率的概念;第二节介绍了P 值的定义和性质;第三节和第四节分别介绍了检验统计量在独立情形和相依情形下FDR 控制的检验方法,第五节介绍了两种真实零假设0m 或比值0π的估计方法,为后面参数混合模型的估计方法奠定了基础。
第3章 参数混合模型和非参数混合模型的估计
3.1 引言
在实际多重假设检验的研究中,我们往往使用随机的检验。
当0=i H 时,统计量X 的密度函数记为()⋅0f , 当1=i H 时,统计量X 的密度函数与某个未知的参数δ有关,记为()δ⋅1f . 这里的m i ,,1⋯=. 如果δ固定,统计量X 的密度函数就可以表示为
()()()(),11000δππx f x f x f -+= (3-1) 与之相对应的P 值密度函数就可以表示为
()()()()().1,0,11000∈-+=p p g p g p g θππ (3-2) 其中上式中的()⋅0g 和()θ⋅1g 分别表示P 值在零假设和备择假设下的密度函数。
显然模型(3-1)和(3-2)是关于δ的参数混合模型。
在模型(3-1)中,参数0π与δ均是可辨别的,其中δ表示冗余参数。
同理,在模型(3-2)里面,参数0π和θ也是可辨别的。
下面分别研究在正态混合分布模型和Beta 混合分布模型下0π的估计方法。
3.2 正态分布混合模型
为了方便研究,本节我们对模型(3-1)中的密度函数加以条件限制。
假设统计量X 在零假设下服从标准正态分布,即()1,0~N X , 那么()⋅0f 为标准正态分布密度函数,我们把它记为()⋅φ; 在备择假设下,统计量()1,~δN X , 也就是说,()δ⋅1f 是期望为δ,方差为1的正态分布密度函数,记为()δφ-⋅。
从而随机变量X 的其密度函数可以表示为
()()()().100δφπφπ--+=x x x f (3-3) 在这个模型中,参数()δπ,0是可辨别的,其中0π是我们要研究的参数,δ为冗余
参数。
对于任何一个样本m x x x ,⋯,,21,如果样本容量足够,就可以由样本的前两阶矩得到方程组
()()()
⎪⎪⎩⎪⎪⎨⎧+-+=-=∑∑==1112001201δππδπm
x m x m i i
m i i
解这个方程组,得 ()
(
)
.1,122
10112
∑∑∑
∑==∧
==∧
--
=-=
m i i
m i i m
i i
m
i i
m
x m x x
m x
πδ (3-4)
即为参数()δπ,0的矩估计。
我们利用基于P 值的最小二乘估计来研究参数()δπ,0的估计方法,这里我们只考虑右侧检验。
令()⋅Φ表示标准正态分布的分布函数,则有()x p Φ-=1,即
()X P Φ-=1。
于是有
()()()ααμμαΦ==<==>0Pr 0Pr H X H P , ()()()δμμααα-Φ==<==>1Pr 1Pr H X H P ,
其中αμ表示正态分布的上侧α分位点,α为检验水平,有()αμα-=Φ1。
记
()}{αα≤=i P R #,则有
()}{αα≤=P m ER Pr
()()()}{1Pr 10Pr 00=≤-+=≤=H P H P m απαπ ()[]()()[]}{1Pr 110Pr 100=>--+=>-=H P H P m απαπ ()()[]}{δμπαπα-Φ--+=1100m
()()}{αμδπαπ-Φ-+=001m 上式可以写成
()()()[].0ααμδαπμδα-Φ-=-Φ-m ER (3-5) 我们取()m i p i ,,2,1,⋯==α,则上式变为
()()()()()()()[]
.0i i p i p i μδΦp m p ER --=-Φ-πμδ
再令
()()()
()()
i i p i p i i m
i
z p y μδμδ-Φ-=
-Φ-=, (3-6)
若δ已知,对点列()m
i i i z y 1,=作最小二乘估计,可以得到参数0π的估计值,即由
2
10min ∑=∧
⎪⎭
⎫
⎝⎛-m
i i i y z π得到参数0π的估计值
.1
2
10
∑
∑==∧=
m i i
m
i i i LS
y z y π (3-7)
而实际上δ是未知的,而可以采用矩估计的方法得到它的初始估计值。
那么这个算法的步骤如下: 算法一:
第一步:采用矩估计方法,由(3-4)式得到参数δ和0π的估计,即参数δ和0π的初值,记为0∧δ和0
0∧π;
第二步:令0
∧=δδ,带入到(3-6)式中,计算点列()m i z y i i ,,2,1,,⋯=;
第三步:对点列()m i z y i i ,,2,1,,⋯=作最小二乘估计,由(3-7)式得到0π新的估计值∧
0π;
第四步:利用2
10∑=∧
⎪⎭
⎫
⎝⎛-m
i i i y z π的最小化方法,求得∧δ;
第五步:令∧
∧=00
0ππ,∧
∧=δδ0
,重复计算第二步至第四步,直到估计值收敛为止。
我们再来利用统计量的拟合方法来研究参数()δπ,0的估计值。
定义
()}{()()()()()
()
,#:λλλλμλλλR X X N R X i N N i i
i ∑∈==≥=,,
并记()()m R M λλ=, 其中λ为给定的检验水平,经过计算得
()()()()()
,1,,0,λλλλλμδμδφδμλμφμ-Φ-+==≥=
=≥H X X E H X X E
()()()
()
λλλμμμ≥=≥⋅==≥=X H X H X H Pr 0Pr 0Pr 0Pr
()()
()λμδπλπλ
πλh ≡-Φ-+=
0001
从而()().11Pr λμλh X H -=≥=
()()H X X E E X X E X H ,λμλμμλ≥=≥≥
()()()()()1,10,=≥-+=≥=H X X E h H X X E h λλμλμλ
()()()[]()()⎥⎦⎤
⎢⎣
⎡-Φ-+⋅-+⋅=λλλμδμδφδλλμφλh h 1
()()()()[]()()
λλλλμδπλπμδφμδδπμφπ-Φ-+-+-Φ-+=
000011 (3-8)
()()()()⎪⎩⎪⎨⎧⎪⎭⎪
⎬⎫⎥⎥⎥⎦⎤⎢⎢⎢⎣
⎡==∑≥λλλλλμr R R X E E X E i X i i :
()
()()⎩⎨⎧⎭⎬⎫
≥⋅=λμλλX X E r r E 1
()λμ≥=X X E (3-9)
()}{()()()λλλμδπλπμμλ-Φ-+=≥=≥=001Pr :#X m
X i E
EM i
那么,由(3-8)和(3-9)式可以得到
()()()()[]()λλμμδπλπλλ≥⋅-Φ-+=⋅X X E X E EM 001
()()()()[]λλλμδφμδδπμφπ-+-Φ-+=001 ()()()λλμϕπμφπ001-+= (3-10)其中()()()λλλμδδμδφμϕ-Φ+-≡。
从而 (3-10)式可化为
()()()λμϕλλ-⋅X E EM ()()[]λλμϕμφπ-=0
与前文类似,令()m i p i ,,2,1,⋯==λ,则有
()()()()()()
.i
i
i
p p p μδδμδφμϕ-Φ+-=
记
()()()
()
()()
()i i i p m
i m j i p p i j x m
z y μϕμϕμφ-=
-=∑+-=1
1, (3-11)
于是,由()2
10∑=-m
i i i y z π可以得到0π的最小二乘估计值同式(3-7) 。
从而这个算法的具体步骤如下: 算法二:
第一步:由(3-4)式得到参数δ和0π的初值0
∧δ和0
0∧π;
第二步:令0
∧=δδ,代入到(3-11) 式中,计算点列()m i z y i i ,,2,1,,⋯=; 第三步:对点列()m i z y i i ,,2,1,,⋯=作最小二乘估计,由(3-7) 式得到0π的新估计值∧
0π;
第四步:利用2
10∑=∧
⎪⎭
⎫
⎝⎛-m
i i i y z π的最小化求得∧δ;
第五步:令∧
∧=00
0ππ,∧
∧=δδ0
,重复计算第二步到第四步到估计值收敛为止。
3.3 Beta 分布混合模型
这一节我们来研究关于P 值的模型(3-2). 由第二章P 值的性质,我们可以考虑采用Beta 分布来拟合模型,那么关于P 值的模型(3-2)转化为
()()()[],1,0,,1,,*
1000∈-+=p b a p g b a p g πππ (3-12)
其中()b a p g ,*
1是参数为b a ,的Beta 分布的密度函数,其具体表示如下:
()()()()
()[].1,0,1,1
1*
1
∈-ΓΓ+Γ=--p p p b a b a b a p g b a
特别情况下,当1=b 时,()(]1,0,1,1*
1∈=-p ap a p g a ,模型(3-12)就转化为
()()(]1,0,1,1000∈-+=-p ap a p g a πππ (3-13) 其中()1,0∈a 。
再来看模型(3-12),我们令()b a G ,*1⋅表示为Beta 分布的分布函数,则有
()().,11Pr *1b a G H P αα-==>
与前面所描述的正态混合模型类似,同样可以采用最小二乘估计。
设α为检验水平,记()}{αα≤=i P R #,则有
()()()[]
.,1*100b a G m ER απαπα-+=
令()m i p i ,,2,1,⋯==α,代入上式中,经过计算得到
()()()()()()()[]
.,,1
*10*1b a p G p b a p G p ER m
i i i i -=-π 记
()()()()(),,,,*1*1b a p G m
i
z b a p G p y i i i i i -=
-= 则若参数b a ,已知,利用点()m
i i i z y 1,=的最小二乘估计方法,可以求得0π的估计值同式(3-7)。
和上节相同,我们仍然采用矩估计方法求得参数的初值。
由样本的前三阶矩
可以得到下列方程组
()()()()()()()()()()()⎪⎪⎪
⎩⎪⎪⎪⎨⎧-++++++++=-+++++=-++=∑∑∑===m i i m i i m i i b a b a b a a a a p m
b a b a a a p m b a a
p m
10031002100
1212141111131
11211ππππππ (3-14)
解这个方程组,得到0,,πb a 的初始矩估计,记为0
000,,∧∧∧πb a 。
我们来研究参数b a ,的极大似然估计方法。
由模型(3-12),其对数似然函数为
()()}
{∑=+=m
i i b a p g b a l 1*
1
10,ln ,ππ, 上式中的011ππ-=. ()b a p g ,*
1关于a 的偏导数为
()()()()()()b a p g p a a b a b a a b a p g i i i ,ln ,*1''*
1
⎥⎦
⎤⎢
⎣⎡+ΓΓ-+Γ+Γ=∂∂ ()()[]()b a p g p a b a i i ,ln *
1+ψ-+ψ=
这里()⋅ψ表示Digmma 函数,即()()()
x x x ΓΓ=ψ',()x Γ为Gamma 函数。
从而有
()()
()∑=∂∂+=∂∂m i i i a b a p g b a p g a b a l 1*
1*
1101
,,,πππ ()()
()()[]∑=+ψ-+ψ+=m
i i i i p a b a b a p g b a p g 1*
110*11ln ,,πππ ()()()[]∑=+ψ-+ψ=m
i i i p a b a p 1
1ln π
其中
()()()()
.,,,,,1Pr *
110*1101b a p g b a p g p b a H p i i i i i πππππ+=== 同理,有
()()()()()[].1ln ,1
1∑=-+ψ-+ψ=∂∂m
i i i p b b a p b b a l π 令
()()0,,0,=∂∂=∂∂b
b a l a b a l ,。