勾建伟, 夏业茂
(南京林业大学理学院, 江苏 南京 210037)
分类数据在社会调查中普遍存在.在客服满意度调查中, 消费者对服务产品的满意度通常存在如下评价:“很满意”,“满意”, “一般”,“不满意”,“非常不满意”; 人口普查问卷通常设计为诸如性别、民族、受教育程度、职业等若干选项用以了解人口和住户的基本情况; 在产品质量检测时对产品质量等级的认定等, 这都形成分类数据.对分类数据的统计分析, 一般着眼于研究外部因素如何影响或作用于“类”的机制, 这通常存在两种办法: 一是基于类的观测频率并结合各种logit模型来进行统计推断.这方面存在大量的文献和著作, 经典的《Categorical data analysis》[1]一书对此做了详细的分析和回顾.但这种方法存在一些问题, 譬如当分类变量个数较多时, 类别总数会急剧增加, 这必将大大提高计算的强度和分析的复杂度, 另外, 该方法在刻画多重分类变量的关联性时不尽如意, 因为它需要借助于偏相关这一数学工具, 这无疑割裂了多重变量的整体关联性.另一种是引入潜变量方法.该方法假定每个分类变量都联系一个连续潜变量或向量, 类的界定是通过潜变量或向量的实现值落在不同区域或窗口来完成的.例如产品的使用期达到一定期限时可以认为是次品, 考试成绩达到一定分数时视为合格等.相比于logit模型而言, 潜变量方法不仅在统计建模和计算方面具有较大优势, 而且在刻画不同分类变量的关联性方面也十分自然: 连续变量之间的相关性在一定程度上就能够体现出分类变量之间的关联性.[2-7]
基于潜变量方法的分类数据分析也存在一些问题, 例如过于依赖潜变量的分布指定等.但最突出的问题是模型或似然对参数存在不可识别性, 也称为模型的不确定性.若一个参数统计模型p(Y,θ)有p(Y,θ1) =p(Y,θ2)蕴含θ1=θ2.则称p(Y,θ)基于Y对未知参数θ是可识别的.[8]一个对参数不能识别的模型将对估计构成极大威胁, 它会产生不相合估计.带有分类变量的潜变量模型的不确定性归根结底是由潜变量的不可观测性所造成的.为了解决分类数据潜变量模型的识别性问题, 很多作者对此作了研究, 但主要工作是建立在单个分类变量的基:上并局限于回归分析领域[9-13].这对于带有多重分类变量的潜变量模型的可识别性显然是不够的.近年来, 部分学者在多元分析领域内对多重二分、次序变量展开统计建模, 并针对模型的可识别性提出了一些具体方法和建议[14-16], 但这些工作大多是面向具体模型且没有从理论上给予严格保证, 很多方法只是保证模型是局部可识别的.
本文主要是对带有多重二分、次序或名义变量的潜变量模型的可识别性在因子分析框架内给出了若干使用方便的充分条件.这些条件本质上是将因子分析模型的可识别性和分类数据的潜变量模型的可识别性结合起来, 但并不仅仅是简单组合, 而是力求从模型结构本身去提出条件.这些条件为模型的统计计算提供了一定的便利.
设对i= 1,··· ,n,yi= (yi1,··· ,yip)T为p维独立的二分观测向量, 其中yij为取值0,1的随机变量,j=1,··· ,p,p ≥2,“T”表示转置.在潜变量模型框架内, 通常假定存在一个p维连续潜变量使得
其中,
μ为p×1维均值向量,Σ为p×p阶协方差矩阵.对应的似然函数为
其中,φ(·|μ,Σ)为Np(μ,Σ)的概率密度函数,
为由yi确定的积分区域.
模型(2.1)-(2.2)通常称为多元Probit模型.然而似然(2.3)却对μ,Σ具有不可识别性.事实上,
定理2.1若存在{μ(1),Σ(1)},{μ(2),Σ(2)}使得
当且仅当存在一个正定对角矩阵D=diag{d1,··· ,dp}(dj >0), 使得
证必要性: 首先注意到样本量的增加并不能改变模型的确定性, 因此考虑n=1, 略去下标i.令由
可知, 若(2.4)成立, 则必然有
其次, 注意到:
其中,R为Σ的相关系数矩阵.上述方程组共计有2p个方程, 由(2.7)及上述方程组蕴含R(1)=R(2), 从而Σ(2)=DΣ(1)D.
充分性: 设(2.5)成立.因为
证毕.
为了更仔细地说明方程组(2.8)蕴含R的唯一性, 我们以p=2为例.不难计算,
其中,φ(·)为一元标准正态概率密度函数, 积分区域
定理2.1表明似然函数Lb并不能由μ,Σ唯一确定, 而是在参数的尺度变换(等价类)意义下唯一确定.在实际应用中, 为了破坏这种同变性, 通常采用参数约束的办法, 这导致如下的结果.
推论2.1若限制μ为分量不为0的已知向量μ0或Σ为相关系数矩阵, 则似然(2.3)关于未知参数是可识别的.
证由假设可知, (2.5)中D=Ip.证毕.
在因子分析领域, 通常将上述的延伸至带有因子变量的因子模型[17]:
其中,Λ为p×m阶因子负荷矩阵,ωi ~Nm(0,Φ)为m×1维因子变量,δi ~Np(0,Ψδ)为p×1维误差向量, 且ωi与δi两者独立.记θ={μ,Λ,Φ,Ψδ}, 则(μ,Σ), 其中,
此时似然函数依然为(2.3)形式, 只不过需要将其中的Σ换成Σ(θ).很明显, 该似然函数关于θ不可识别, 其原因有两种: 1)因子模型本身不确定; 事实上注意到(2.11)中的边际分布在因子可逆变换下会保持不变.2) Probit模型(2.1)不确定性; 这由定理2.1给出.为了保持最终的似然关于未知参数是可识别的, 一种有效的方法是在保证因子模型确定的基:上利用定理2.1通过约束参数来保证最终似然对参数具有识别性.
对因子分析模型的模型确定性研究已经有相当长的历史, Bollen[17]等就对此就作出过深刻的研究.但到目前为止, 仍无一个使用方便的充分必要条件来确定该模型.大多数的做法都是面向具体模型通过约束参数来达到模型的确定性.本文引用Bollen[17]的两指标准则.该结果虽然是一个充分条件, 但它具有较大的普适性: 对一般的因子模型(无论正交或倾斜因子)均适合, 而且因为它对因子负荷矩阵和误差协方差矩阵的元素分开进行约束, 这在计算和理论分析上具有较大优势.我们引用如下:
两指标准则[17]:
1) 每个因子都尺度化, 即对每个因子ωik, 存在λjk=1; 亦即Λ的每一列都存在一个元素1;
2)Λ的每一行都有一个非零元素;
3) 联系每个因子的指标个数至少为2; 即Λ的每一列至少存在两个非零元素;
4)Ψδ为对角正定矩阵;
5)Φ正定, 且至少存在一对非0的非对角元Φjk0().
对上述条件的说明是必要的: 条件1)是对因子进行标度.由于因子变量为不可观测变量,其单位一般不明确, 限定因子系数为1是将因子变量的尺度与观测变量的尺度等同起来, 这便于解释该因子对其它指标变量的贡献; 条件2)则表明每个指标变量都联系一个因子, 这也是必要的, 因为指标的主要作用是明示因子; 条件3)要求明示因子变量时需要多重而非单个指标,这主要是因为因子本身就是来解释多指标的相关性; 条件4)则要求因子变量能够解释指标间的全部相关性; 条件5)要求至少存在一对因子不独立.这个条件可以适当放宽.如果要求“每一列非零元素的个数至少为3”时, 条件5)的因子可以允许相互独立.
下面我们考虑模型(2.1), (2.11)的可识别性问题.根据推论2.1, 在保证因子模型确定的基:上, 只需要约束Σ(θ)为相关系数矩阵即可.但这种约束破坏了因子负荷参数与协方差参数的分离性, 将会给统计计算造成一定的困难.下面的结果给出了一个充分条件.该条件只要求在模型(2.11)确定的基:上约束误差的协方差阵为单位矩阵即可, 这保持了参数的分离性, 因此在计算上较为方便.
定理2.2假定模型(2.11)中的Λ,Ψδ和Φ满足两指标准则, 且Ψδ=Ip, 则边际似然(2.3)关于μ,Λ,Φ是确定的.
证由定理2.1可知, 若Lb(θ(1)) =Lb(θ(2)), 则存在正定对角矩阵D= diag{d1,··· ,dp}(dj >0), 使得
由Σ的确定性立即可知d1=···=d4=1.0.
当然,在(2.11)假定下,似然函数的不确定性也可以通过限制均值μ=μ0来完成,其中μ0为分量不为0的已知向量, 但这需要事先明确或估计出μ0的位置, 这又要归结到Σ约束上.
若定理2.2中的Φ含有参数结构Φ(φ) 且Φ(φ)关于φ确定时, 则最终的似然函数也关于φ确定.
推论2.2假定模型(2.11)中的Λ,Ψδ和Φ满足两指标准则, 且Φ=Φ(φ).若限定Ψδ=Ip且Φ(φ)关于φ确定, 则边际似然(2.3)关于μ,Λ和φ确定.
推论2.2一个直接应用就是结构方程模型的确定.在潜变量分析领域, 结构方程模型[17]是一种用来描述多重因子之间内部相互关系和因果关联的统计方法.该模型将因子变量ωi分为m1个内在因子ηi和m2= (m-m1) 个外在因子ξi.内在因子之间的关联性以及外在因子对内在因子的影响方式是通过如下结构方程来体现:
其中, 误差变量ζi与外在因子ξi以及模型(2.11)中的误差变量δi相互独立,B为主对角元素为0的m1×m1阶矩阵,Γ为m1×m2(m2=m-m1)阶回归系数矩阵,Ψζ为m1×m1阶对角正定矩阵,Ξ >0.
推论2.3假定模型(2.11)中的ωi={ηi,ξi}满足(2.15)-(2.16).若Λ满足二指标准则,Ψδ=Ip,B为下三角均阵, 则似然(2.3)完全由参数μ,Λ,B,Γ,Ψζ和Ξ确定.
证很显然, 在(2.15)-(2.16)假定下,ωi ~N(0,Φ), 其中,
B0=Im1-B.不难验证, 若B为下三角阵(注意其对角线元素为0), 则Φ关于B,Γ,Ψζ和Ξ确定.由推论2.1, 结论立得.证毕.
关于Φ确定性的一个具体例子为
此时,
显然,Φ关于ξ是确定的; 其次注意到
由此可知Φ对Γ11,ψζ1,β21,Γ21,ψζ2确定, 从而对B,Γ,Ψζ,Ξ确定.证毕.
注意, 当推论2.3中B=0时结论自然成立, 这是无递推结构方程模型.
将多重二分变量推广为多重次序变量, 则得到带有多元次序变量的潜变量模型.为简单起见, 设yi= (yi1,··· ,yip)为p维观测向量, 其中yij均为取值于{0,1,··· ,K}(K ≥1)的次序变量.K=1即为二分变量.需要指出是, 式(2.1)式中0为的门限值并非必要.事实上, 0可以用本节的αj替代.为方便起见, 本节假定K ≥2.在潜变量模型框架内, 假定yi联系一个p维的连续潜变量(μ,Σ), 使得
其中,-∞=αj0<αj1<···<αjK <αjK+1=+∞为用来界定类的门限值.记
则似然函数为
其中,
是由yi确定的积分区域.
与二分变量类似, 似然(3.3)对α,μ,Σ存在不可识别性.
定理3.1若存在{α(1),μ(1),Σ(1)},{α(2),μ(2),Σ(2)}满足
当且仅当存在一个对角正定矩阵D=diag{d1,··· ,dp}(dj >0), 使得
其中αℓ为α的第ℓ列.
证类似地, 我们考虑n=1.令首先考虑必要性.注意到yi的边际分布为
其中,Φ(·)为标准正态分布函数.由(3.4)可知, 对ℓ=1,··· ,K,j=1,··· ,p,
其次, 注意到
其中R为Σ的相关系数矩阵.(3.4),(3.7)及上式蕴含R(1)=R(2).
令T(j)=则T(2)=DT(1).又因为
从而Σ(2)=DΣ(1)DT, 因此式(3.5)成立.
反之, 若(3.5)成立, 则在(3.4)的左边积分中令(1)=D-1((2))立得右式.证毕.
定理3.1表明似然函数Lp是在参数的仿射变换(等价类)意义下唯一确定.同样地, 在实际应用中, 为了破坏这种同变性, 通常采用约束参数的办法, 这导致如下结果.
定理3.2当(μ,Σ)时, 若下列条件之一成立:
1) 若对j=1,··· ,p, 约束αj1,αjK为固定值;
2) 若对某个ℓ, 约束αℓ为固定值, 且Σ为相关系数矩阵;
3) 若约束μ为固定值,Σ为相关系数矩阵;
4) 若对某个ℓ, 约束αℓ为固定值, 且μ为固定值.
则边际似然(3.3)关于自由未知参数是可识别的.
证由定理3.1可知,若(3.4)成立,则存在一个对角正定矩阵D=diag{d1,··· ,dp}(dj >0),使得
1) 由假设可知
由上述定理可知, 为了保证似然Lp对未知参数的可识别性, 至少需要在门限值、均值、协方差矩阵三组参数中约束两组参数.在回归分析领域中, 文献上常采用定理3.2中的约束1)或2)[14-15]; 这就需要事先知道门限的具体位置.一个可行的做法是: 令αj1=Φ-1αjK=Φ-1(, 其中fj1,fjK分别为yij <1 和yij <K的频率,Φ-1(·)为标准正态分布函数的反函数.
类似地, 我们将模型(3.1)中的推广为因子模型(2.11), 此时Σ=ΛΦΛT+Φ.与前面二分变量类似, 此时边际似然存在两方面不确定性.我们依然假定Λ,Ψδ,Φ满足两指标准则.
推论3.1假定模型(2.11)中的Λ,Ψδ和Φ满足两指标准则, 且下列条件之一成立:
1) 若对j=1,··· ,p, 约束αj1,αjK为固定值;
2) 若对某个ℓ, 约束αℓ为固定值, 且Ψδ=Ip;
3) 若约束μ为固定值, 且Ψδ=Ip;
4) 若对某个ℓ, 约束αℓ为固定值, 且μ为固定值.则边际似然(3.3)关于α,μ,Λ,Ψδ,Φ是可识别的.
证首先, 两指标准则保证因子模型(2.11)是可识别的.
1) 由定理3.2的1), 立得结论;
2) 由定理3.1可知,Ψδ=Ip蕴含(3.5)中的D=Ip.再由定理3.2中2)可知似然函数关于均值参数和门限值是可识别的;
3) 同2);
4) 由定理3.2的4), 结论成立.
证毕.
类似地, 若ωi满足结构方程(2.15)-(2.16), 则有下列结果:
推论3.2假定模型(2.11), (2.15)中的Λ,Ψδ和Φ满足两指标准则,B为下三角阵, 且下列条件之一:
1) 若对j=1,··· ,p, 约束αj1,αjK为固定值;
2) 若对某个ℓ, 约束αℓ为固定值, 且Ψδ=Ip;
3) 若约束μ为固定值, 且Ψδ=Ip;
4) 若对某个ℓ, 约束αℓ为固定值, 且μ为固定值.
则边际似然(3.3)关于α,μ,Λ,Ψδ,B,Γ,Ψζ,Φ是可识别的.
证结合推论2.3和定理3.2, 不难得出上述结论.证毕.
名义变量又称为无次序变量.与次序变量类似, 名义变量取值依然表现为多岐性, 但此时yij ∈{0,1,··· ,K}的取值仅仅代表类别, 而无次序含义.譬如, 在若干种商品组成的选项中, 很难界定一种商品比另一种商品优或劣, 这时就表现出名义属性来.此时用潜变量模型(3.1)来界定这样的类并不恰当, 但可以借助于该思想, 对每个yij, 引入一个K维连续潜在向量使得
其中,
μ为pK×1维截距向量,Σ >0为pK×pK阶正定协方差矩阵.在某些情形下, 为了降低参数维数, 可以取Σ=U ⊗V, 其中U,V分别为p×p和K ×K阶正定矩阵, 这导致服从矩阵正态分布.
此时似然(4.3)完全由参数{μjk/σjk,σjk/σjℓ,Rjk}或等价地由参数{μjk/σj1,σjk/σj1,Rjk}确定.因此, (4.4)蕴含着
与多元次序变量潜变量模型识别性类似, 可以考虑限制参数来保证似然对参数的可识别性.
推论4.2对于模型(4.1)和(4.2), 记
若限制第ℓ(= 1,··· ,K)列μ(ℓ)元素为不为0的固定值或对j= 1,··· ,p, 限制Σjj的某个对角元素为1, 或trΣjj=dj0>0为常数, 则似然(4.3)关于未知自由参数是确定的.
我们首先考虑模型(4.1), (4.11)似然的可识别性.
定理4.2考虑因子模型(4.11).若Λj,Ψδj,Φ满足两指标准则, 且对j=1,··· ,p,ψδj,11=1, 则似然(4.3)对参数μj,Λj,Ψδj,Φ可识别.
证由定理4.1可知,若(4.4)成立,则存在一个对角正定矩阵D=diag{d1,··· ,dp}(dj >0),使得
特别地,
类似于定理2.2证明, 可以调整Λ(1)和Φ(1)使得
证毕.
定理4.2要求(4.11)中Λj,Φ,Ψδj都满足二指标准则, 这在实际应用中往往难以实现.例如m=1, 至少K=2 当m=2时,K至少要等于3.可以将上述条件稍微放宽一点, 得到如下结果.
定理4.3考虑因子模型(4.11).记
若Λ,Φ,Ψδ满足二指标准则, 且对j=1,··· ,p,ψδj,11=1, 则似然(4.3)对参数μj,Λj,Ψδj,Φ可识别.
证由定理4.1可知, 存在D使得
类似与定理2.2证法, 我们有
一个符合定理4.3但不符合定理4.2的例子为: 考虑p=K=m=2,
显然,Λj并不符合两指标准则, 但Λ却是符合的.
不难验证,Σ关于λjk,φjk和ψδj是确定的.
我们按照上述模型产生数据集.总体参数真实值设置为:μ= 0.8×18,λjk= 0.8,ψδj=1.0, (αj1,αj2,αj3)=(-1.6,0.0,1.4)(j=3,4),
对于当前模型而言, 由于共计存在6层2×2×4×4×3×3=576个格子, 为了避免格子中的样本稀疏, 我们取样本容量为6000, 平均到每个格子的频数约为10.
我们对上述模型进行贝叶斯推断.与频率推断相比, 贝叶斯方法不必依赖于大样本理论,且计算简单.考虑如下先验:
为了保证模型的可识别性, 我们将Ψδ取定为单位矩阵.在计算上, 我们执行马尔可夫链蒙特卡洛(MCMC)抽样方法[19].马尔可夫链蒙特卡洛要求从p(μ|···),p(Λ|···),p(α|···) 以及p(Φ|···)中抽样.可以仿照[7], 不难得出这些条件分布的具体形式, 这里从略.我们考虑三组不同的初始值来确定估计的收敛性.收敛性通过EPSR值[18]来判断.图5.1分别给出了各个参数在先验I及不同初始值下EPSR值对迭代次数的关系图.由此图看出, 在1000步以内, 所有的参数的EPSR值都小于1.2.我们截取3000步以后的样本来计算估计偏差、均方误和标准误差.结果由表5.1给出.从表5.1可以看出, 两者估计差别区别并是不是很大, 这说明估计对超参数值的选择具有一定的稳健性, 但Λ12和Φ11估计的均方误和标准误差较大, 这主要是因为二分数据提供的总体信息较少.我们将样本量提高到10000, 则结果明显改善.为节省篇幅, 计算结果不再一一列出.
模型或似然的可识别问题并不仅仅限于本文所讨论的范围, 在诸如回归分析、方差分析、缺失数据分析等领域都存在一定的模型确定性问题.本文主要是针对多重分类数据的潜变量模型的可识别性问题进行了讨论, 给出了若干使用方便的充分条件.这些条件可视为将因子分析模型的可识别性和分类数据的潜变量模型的可识别性结合起来, 但也并不仅仅是简单组合,而是力求从模型结构本身去提出条件.尽管这些只是充分条件而非必要的, 但它们在模型的统计计算方面提供了较大的便利.
图5.1 先验I下各参数EPSR值对迭代次数的关系图: Ψδ =I8
表5.1 两种先验下各参数估计: 样本量N =6000, Ψδ =I8
不可识别的模型并非要被弃用的模型.模型可识别性只是对参数估计构成影响, 并不影响模型对数据拟合评价, 这从似然的不变性不难看出.不可识别的模型虽然对某些参数估计的确构成困难, 但它们的某些函数却是可估的.例如不确定的因子分析模型虽然对Λ,Φ,Ψδ估计造成麻烦, 但并不影响Σ=ΛΦΛT+Ψδ的估计.有些情形下, 不可识别的模型可能有助于快速找到模型解, 如正交因子模型中的因子旋转.特别在计算上, 著名的PXEM算法(PXDA)算法正是利用模型的不确定性来扩大解的搜索范围从而便于快速地进行求解(见文[20]).
最后需要指出的是, 对参数约束只是为了保证非约束参数的估计能够到相合到唯一值, 但这并不意味着得到的估计就一定相合到参数真实值.实际上, 对于不可识别模型而言, 模型只能够识别到可识别的那些参数.例如, 对于正态模型N(α+β,1)而言, 模型只能对α+β识别.为了识别α, 可以将β固定成0, 也可以固定为其它值, 但这两种不同约束导致α的估计相合到的结果一般并不相同.对于那些非约束参数而言, 不同的约束实际上代表不同的模型机制; 估计对约束值也相当敏感, 一般并不具有稳健性.这就引起一个问题: 如何选择一个较为适宜的约束.在本文中, 我们尽管讨论了约束某些参数可以保证模型是确定的, 但没有指明这些约束值是否一定合理, 如何约束需要考虑到问题的实际背景.比如推论2.1中, 将协方差阵约束为相关系数矩阵, 这不并影响潜变量关性的界定.在定理2.2、推论3.1和定理4.2中, 我们将Ψδ约束为一个单位矩阵, 但在实际问题中如果潜变量的协方差阵Σ=ΛΦΛT+Ψδ特征根过小, 上述约束可能导致ΛΦΛT奇异, 其结果是Λ,Φ中非约束参数可能出现不合理的解.一个较好的做法是将Ψδ约束为c0I形式, 其中c0=0.1或0.01.如何对约束参数赋值留作进一步的研究和讨论.