序约束下单向分类方差分析模型的Bayes变量选择

史海芳, 李 聪, 姬永刚

(1. 中国民航大学 理学院, 天津 300300; 2. 吉林大学 数学学院, 长春 130012)

在剂量反应实验中, 时间可作为解释变量, 人们可能在收集数据前就已获得了潜在参数满足序约束的先验信息. 将这些先验信息应用到统计推断中可以提供更有效的估计, 特别是在具有低信噪比的小样本情形下. 文献[1-2]对序约束统计方法进行了全面阐述.

目前, 利用Bayes方法解决约束条件下的统计推断问题已取得许多成果. 文献[3]借助Savage-Dickey密度比方法[4]讨论了序约束下单向分类方差分析(ANOVA)模型的假设检验问题, 利用该方法可近似计算各潜在的假设Bayes因子, 并选择后验概率最大的假设作为最优假设； 文献[5]将该方法推广到其有不等式约束条件的广义线性模型中, 但Savage-Dickey密度比方法需要在每次Markov链Monte Carlo(MCMC)迭代中都近似计算标准正态分布的密度和分布函数值, 计算时间较长； 文献[6-7]通过引进指标变量, 并假设指标变量与参数相互独立, 分别将文献[8]中的方法应用到单向分类和双向分类方差分析模型中, 考虑了序约束下的变量选择问题, 其优点是没有过多的调节参数, 但文献[9]指出若参数先验分布的假设太模糊, 则可能会导致Markov链收敛较慢; 文献[10-11]提出了另一种常用的Bayes变量选择方法SSVS(stochastic search variable selection)方法, 目前SSVS已被扩展到很多模型和应用中, 例如遗传数据分析[12]、 向量自回归模型[13]和分组约束模型[14]等. 该方法的缺点是需要进行调节的参数较多. 为解决上述问题, 文献[15]提出了一种改进的SSVS方法----NMIG(normal mixture of interse Gamma distributions)方法, 并将其应用到基因分析中. 本文将这两种SSVS方法推广到单向分类方差分析模型中， 并考虑序约束下的变量选择问题. 基于MCMC方法引入指标变量, 计算潜在模型的后验概率.

1 再参数化ANOVA模型中的变量选择

令

yij=μi+εij,j=1,2,…,ni,i=1,2,…,k

(1)

表示一个单向ANOVA混合模型, 其中yij表示第j个个体第i次治疗的反映变量,μi表示第i次治疗的平均效应,εij表示误差项且εij～N(0,σ2).实际应用中一般有如下几种序约束：

(i) 简单半序约束μ1≤…≤μk;

(ii) 简单树序约束μ1≥μi,i=2,…,k-1;

(iii) 伞序约束μ1≤…≤μg≥μg+1≥…≥μk, 其中g已知.

对于简单半序约束μ1≤…≤μk, 如果令δm-1=μm-μm-1(2≤m≤k), 则第二个均值μ2可表示为μ1+δ1, 第三个均值可表示为μ1+δ1+δ2,…, 第m个均值(2≤m≤k)可表示为μ1+δ1+…+δm-1,μm的简单序约束可等价地表示为δm的非负约束δm≥0,m=1,2,…,k-1.如果令δ=(δ1,δ2,…,δk-1), 则满足简单半序约束的模型(1)可表示为

y=1nμ1+xδ+ε,δ1≥0, …,δk-1≥0,

(2)

这里1m和0m分别表示元素皆为1和0的m×1维列向量.若均值满足简单树序约束或伞序约束, 则本文也可定义相应的设计矩阵x及再参数化向量参数θ.这些不等式约束可参见文献[3,5].尽管本文给出的方法适用于任何可表示为模型(2)的参数约束, 但为方便, 本文主要考虑简单半序约束.

为能同时进行变量选择和参数估计, 本文考虑两种Bayes变量选择方法： SSVS方法[10-11]和NMIG方法[15]. 本文将SSVS方法扩展到带有序约束条件的单向ANOVA模型中.引入指示变量γi, 并假设γi和δi满足下列先验分布:

δi|γi,φi～TN(0,c(γi)φi,0,+∞),

(3)

(4)

其中:

(5)

π(σ2)∝1,

文献[16]将该先验分布应用于无约束线性模型的变量选择问题中.

上述先验和超先验分布对于其各自参数和超参数均为共轭的, 因此利用Gibbs抽样方法易得参数的后验分布.下面分别给出这两种方法所有参数的满条件分布.

1.1 SSVS方法

1)σ2的条件后验分布.由于

(6)

因此

σ2|y,x,μ1,{γl},{δl}～IGamma(n/2-1,(y-1nμ1-xθ)′(y-1nμ1-xθ)/2).

2)μ1的满条件后验分布.由于

3)γi的满条件后验分布.令γ-i=(γ1,…,γi-1,γi+1,…,γq), 可证明γi的满条件分布服从如下Bernoulli分布:

(7)

其中

ci=f(y|{δi},σ2,γi=1,γ-i)f(δi|γi=1,γ-i)f(γi=1,γ-i),

di=f(y|{δi},σ2,γi=0,γ-i)f(δi|γi=0,γ-i)f(γi=0,γ-i).

4)πi的满条件后验分布.由πi和γi的先验可计算πi的后验分布为

πi|y,x,{βi},σ2,{γi},μ1～Beta(a1+γi,b1-γi+1).

5)δp的满条件后验分布.由于

1.2 NMIG方法

2 数值模拟

下面用数值模拟验证本文方法的有效性, 并与文献[3]和文献[6]提出的方法进行对比, 将这两种对比方法分别简记为Oh方法[3]和Otava方法[6]. 假设k=4, 且均值满足简单半序约束, 则有下列8个备选模型：

(8)

根据文献[3], 假设误差项服从独立的标准正态分布, 考虑以下3种模拟:

模拟1:μ=(0,0,0,0)′; 模拟2:μ=(0,0,0,1)′; 模拟3:μ=(1,2,3,4)′.

设样本个数n=10,30,100, 每次产生10 000个Gibbs样本, 为保证收敛, 丢弃前面的3 000个抽样值.每种固定效应和样本个数的取值都重复500次.表1～表3列出了所有模型的平均后验概率.其中标*处对应真实模型的后验概率.

表1 模拟1中不同方法所得模型的平均后验概率

表2 模拟2中不同方法所得模型的平均后验概率

表3 模拟3中不同方法所得模型的平均后验概率

由表1～表3可见, 4种方法给出的真实模型后验概率均为最大. 一般地, NMIG和Otava比其他两种方法更倾向于提供真实模型更高的平均后验概率. 为评价这些方法在识别正确模型方面的效果, 类似传统假设检验的优势, 本文计算了每种方法选择正确模型的频率P, 结果列于表4～表6. 由表4～表6可见, 在模拟3中Oh方法效果更好, 但在模拟1和模拟2中另外3种方法的正确模型识别率更高. 数值模拟结果表明, 在多数情形下, NMIG方法比SSVS方法效果更好, 这可能是因为在假设参数φi服从Gamma分布的条件下调节参数比固定其为常数更好.

表4 模拟1中不同方法选择正确模型的频率

表6 模拟3中不同方法选择正确模型的频率

3 实例分析

下面用本文方法分析一组由文献[3]分析的实际数据. 实验人员在18个月内测量了20名男孩在8岁、 8.5岁、 9岁、 9.5 岁时的支骨高度(ramus bone heights). 令μ1,μ2,μ3，μ4分别表示4次观测的平均支骨高度. 根据先验信息, Oh方法假设平均支骨高度满足简单半序约束μ1≤μ2≤μ3≤μ4, 并利用Savage-Dickey密度比方法进行了分析. 本文利用SSVS方法和NMIG方法对该数据进行分析, 并与Oh方法和Otava方法进行比较.考虑式(8)模型H0～HF, 假设先验的超参数取值与数值模拟相同, 并产生40 000个Gibbs样本, 将前20 000个样本作为初始值. 表7列出了所有模型的后验概率. 由表7可见, NMIG方法和Otava方法选择了相同的最大后验模型, 同时4种方法均认为8.5岁时的平均支骨高度与9岁时的平均支骨高度有差异.表8和表9分别列出了参数δi(i=1,2,3)的后验均值和可信区间.由表8和表9可见, Otava方法给出的后验均值更大, 而SSVS方法给出了长度更短的Bayes可信区间, 表明SSVS方法在该数据集中表现较好. 图1和图2分别为SSVS方法和NMIG方法中参数δi(i=1,2,3)的后验密度估计.由图1和图2可见, 两种方法都很好地将参数δi(i=1,2,3)控制在非负区间上.

综上所述, 本文讨论了序约束下单向分类方差分析模型中的变量选择问题, 提出了两种基于SSVS方法序约束下的Bayes变量选择方法, 并利用数值模拟和实例分析验证了本文方法的有效性. 本文提出的后验抽样方法较简单, 且易操作.

表7 支骨高度数据中模型的后验概率

表8 支骨高度数据中参数的后验均值

表9 支骨高度数据中参数的95%可信区间

图1 SSVS方法参数δ1,δ2,δ3的后验密度估计Fig.1 A posteriori density estimation of parameters δ1,δ2,δ3 in SSVS method

图2 NMIG方法参数δ1,δ2,δ3的后验密度估计Fig.2 A posteriori density estimation of parameters δ1,δ2,δ3 in NMIG method