统计推断常常通过建立各种模型来探讨因变量与相关协变量之间的关系.关于模型中参数的估计, 例如, 极大似然估计, 最小二乘估计等, 往往涉及到求解一些目标函数的极值问题.以下我们将着重关注极大似然估计的求解问题.很多模型下, 似然方程没有解析解, 这时需要采用数值方法求解参数的估计值, 例如, Newton-Raphson算法, 随机梯度下降法^[1]等.其中, Newton-Raphson算法是最为常用的数值方法之一, 然而, Newton-Raphson算法在实际应用中可能遇到下述问题:每次迭代均需要计算复杂的二次导数矩阵; 当参数维数较大时, 算法涉及到大维矩阵求逆的问题; 由于不同初值的选取, 算法可能收敛到似然函数的极小值点或者鞍点; 多次迭代导致矩阵奇异, 迭代无法继续进行; 由于算法没有上升性, 可能会导致不收敛的情况发生^[2].因此, 发展和研究更为稳定的数值算法是统计研究中的热点问题之一.

近年来, 对Minorization-Maximization算法(以下简称“MM算法”)的研究与应用越来越广泛. MM算法的中心思想是“优化转移”, 关键之处在于构造一个代理函数(surrogate function), 从而将求解复杂的目标函数的极值问题转移为求解具有优良性质的代理函数的极值问题, 进而有效地避免上述提到的Newton-Raphson算法可能遇到的问题. MM算法近来被广泛地应用到统计学各个研究领域, 例如, 稳健回归^[3], 关联度分析^[4], 最小二乘估计^[5-7], 医学影像分析^{[8, 9]}, 分位数回归^[10], 生存分析^[11], DNA排列分析^[12], 多重比较^[13], 变量选择^[14], 判别分析^{[15, 16]}等等. MM算法是一种算法的思想, 而非局限于一种具体的计算方法.算法原理的关键处是构造一个合理的代理函数, 常用方法有:线性化; 参数分离; 规避求解大维矩阵; 利用光滑代理函数替代非光滑目标函数等等^[17-22].

众所周知, 线性回归模型是定量分析研究中最为流行的统计分析方法.然而, 当因变量为分类变量而不是连续变量时, 线性回归模型就不再适用了.对于二元因变量的研究, Logistic回归模型是使用最为普遍和广泛的分析方法.在Logistic回归模型下, Böhning & Linday (1988) 为其极大似然估计的求解提出了一种二次下界(Quadratic Lower-Bound)算法(以下简称“QLB算法”)来构造MM算法中的代理函数^[2]. QLB算法用一个与参数无关的二次下界矩阵替代了负的观测信息矩阵, 从而, 在求解极大似然估计的整个迭代过程中仅需要计算一次该二次下界矩阵的逆, 而非如Newton-Raphson算法在每次迭代中均需要计算一次观测信息矩阵的逆.

本文以探讨Logistic回归模型中参数极大似然估计问题的QLB算法理论为支撑, 主要研究其在实际中的应用.首先, 我们利用数值模拟方法比较了Newton-Raphson算法和QLB算法, 并评估了极大似然估计的大样本性质在有限样本下的表现.进一步, 利用QLB算法分析了三个实际数据:急性淋巴细胞性白血病病人数据, 教职工升职情况数据, 农村居民健康行为数据.本文结构如下:第2节中, 我们首先介绍了一般情形下的QLB算法的原理和方法, 接着, 探讨在Logistic回归模型下, 如何利用QLB算法求解参数极大似然估计.第3节中, 我们先利用数值模拟比较Newton-Raphson和QLB两种算法, 最后应用QLB算法分析了三个实际数据.

本节先介绍一般情形下的QLB算法的原理和方法, 然后, 在Logistic回归模型下, 利用QLB算法求解回归参数的极大似然估计.

设 $l(\theta)$为感兴趣的对数似然函数, $\theta$为待估的 $q$维参数, 则 $\theta$的极大似然估计为

$ \begin{eqnarray*} \label{eq1} \hat{\theta}=\mathop{\text{argmax}}_{\theta\in \Theta} l(\theta). \end{eqnarray*} $

(2.1)

上述极大值问题的求解常常转化为求解得分方程:

$ \frac{\partial l(\theta)}{\partial \theta}=0. $

然而, 实际中, 上述方程的解通常没有解析表达式, 故而需要采用数值方法求解方程, 其中Newton-Raphson算法是广泛使用的方法之一.设 $\theta^{(n)}$表示第 $n$步迭代所求得的 $\theta$的解, Newton-Raphson迭代公式为

$ \begin{equation} \theta^{(n+1)} = \theta^{(n)} - \left(\frac{\partial^2 l(\theta^{(n)})}{\partial \theta \partial \theta'}\right)^{-1} \frac{\partial l(\theta^{(n)})}{\partial \theta}. \end{equation} $

(2.2)

前面提到, 在实际应用中, Newton-Raphson算法常会遭遇初值选取, 观测信息阵奇异以及算法可能不收敛等问题.而此时, MM算法可以有效地避免上述问题的发生, 故而成为一种运用越来越广泛的方法.

MM算法的主要思想是“优化转移”, 其关键之处在于构造一个代理函数 $Q(\theta|\theta^{(n)})$, 使之满足

$ l(\theta)-Q(\theta|\theta^{(n)})\geq 0, \qquad\qquad \forall\theta\in\Theta, $

且等号成立当且仅当 $\theta=\theta^{(n)}$.若定义

$\begin{eqnarray*} \theta^{(n+1)}=\mathop{\text{argmax}}_{\theta\in \Theta} Q(\theta|\theta^{(n)}), \end{eqnarray*} $

则可得如下事实:

$ \begin{eqnarray*} l(\theta^{(n+1)}) &=& l(\theta^{(n+1)})-Q(\theta^{(n+1)}|\theta^{(n)})+Q(\theta^{(n+1)}|\theta^{(n)}) \\ &\ge&l(\theta^{(n)})-Q(\theta^{(n)}|\theta^{(n)})+ Q(\theta^{(n)}|\theta^{(n)})\\ &=& l(\theta^{(n)}). \end{eqnarray*} $

此性质使得MM算法与EM算法一样具有上升性, 这一优良性质.正是由于这种上升性, 使得(2.1) 中的极值问题可以转化为

$ \begin{eqnarray*}\label{eq3} \hat{\theta}=\mathop{\text{argmax}}_{\theta\in \Theta} Q(\theta|\theta^{(n)}). \end{eqnarray*} $

(2.3)

即将求解似然函数 $l(\theta)$极大值点的问题转移为求解代理函数 $Q(\theta|\theta^{(n)})$极大值点的问题.

构造代理函数的方法很多, 文献[17-22]以及文献[11]等分别对不同的统计模型构造了与之相应的代理函数. Böhning & Linday (1988) 提出了一种QLB算法来构造代理函数^[2].假设 $l(\theta)$满足二阶连续可微, 将 $l(\theta)$在 $\theta^{(n)}$处进行Taylor展开, 可得

$ l(\theta) = l(\theta^{(n)}) + \left(\frac{\partial l(\theta^{(n)})}{\partial \theta}\right)' (\theta-\theta^{(n)}) + \frac{1}{2}(\theta-\theta^{(n)})' \left(\frac{\partial^2 l(\tilde\theta)}{\partial \theta \partial \theta'}\right) (\theta-\theta^{(n)}), $

其中 $\tilde\theta$为 $\theta$和 $\theta^{(n)}$连线上的一点. QLB算法关键之处在于找到一个与 $\theta$无关的负定矩阵 $B$, 使得

$ \frac{\partial^2 l(\theta)}{\partial \theta \partial \theta'} \geq B, \qquad\qquad \forall \theta\in\Theta. $

此时, 可构造代理函数形式如下:

$ Q(\theta|\theta^{(n)}) = l(\theta^{(n)}) + \left(\frac{\partial l(\theta^{(n)})}{\partial \theta}\right)' (\theta-\theta^{(n)}) + \frac{1}{2}(\theta-\theta^{(n)})'B(\theta-\theta^{(n)}). $

注意到, 上述代理函数 $Q(\theta|\theta^{(n)})$用 $B$替代了似然函数 $l(\theta)$二阶展开式中的 $\left(\frac{\partial^2 l(\tilde\theta)}{\partial \theta \partial \theta'}\right)$, 且满足

$ l(\theta)-Q(\theta|\theta^{(n)})\geq 0, \qquad\qquad \forall \theta\in\Theta, $

且 $l(\theta^{(n)})=Q(\theta^{(n)}|\theta^{(n)})$.故而可构造求解极值问题(2.3) 的迭代公式如下:

$ \begin{equation} \theta^{(n+1)} = \theta^{(n)} - B^{-1} \left(\frac{\partial Q(\theta^{(n)}|\theta^{(n)})}{\partial \theta}\right) = \theta^{(n)} - B^{-1} \left(\frac{\partial l(\theta^{(n)})}{\partial \theta}\right). \end{equation} $

(2.4)

因为 $\theta^{(n+1)}=\mathop{\text{argmax}}\limits_{\theta\in \Theta} Q(\theta|\theta^{(n)})$, 即每一步的迭代都使得 $Q(\theta|\theta^{(n)})$达到极大值, 故而所得到的估计值序列 $\{\theta^{(n)}:\ n\geq 1\}$会逐渐逼近 $l(\theta)$的极大值点. Böhning & Linday (1988) 给出了上述QLB算法上升性并证明了其算法的收敛性^[2].注意到, 与Newton-Raphson算法相比, QLB算法的最大优点在于:由于二次下界矩阵 $B$与 $\theta$无关, 因此, QLB算法不需要在每次迭代中计算观测信息矩阵的逆, 这个优势在计算大维矩阵的逆时尤为突出.当然, 二次下界矩阵 $B$的选取不是唯一的, 不同的 $B$会产生不同的迭代步长, 从而产生不同的收敛速度.

Logistic回归模型广泛地应用于各个领域, 例如, 流行病学, 生物医学, 经济学等等.假设有 $N$个独立样本, 记 $Y_i$为第 $i$个样本的二元响应变量, $X_i$为第 $i$个样本的 $q$维协变量, $i=1, \cdots, N$.我们考虑下述Logistic回归模型:

$ P(Y_i=1)=\frac{e^{X_i'\theta}}{1+e^{X_i'\theta}}, \quad i=1, \cdots, N, $

其中 $\theta$为待估参数.记 $p_i=P(Y_i=1)$, 对数似然函数有如下形式:

$ \begin{eqnarray*} l(\theta) &=& \sum\limits_{i=1}^{N}\left(Y_{i}\log p_{i}+(1-Y_{i})\log(1-p_{i})\right)\\ &=&\sum\limits_{i=1}^{N}\left(Y_{i}(X_{i}'\theta)-\log(1+e^{X_{i}'\theta})\right). \end{eqnarray*} $

从而得分方程和观测信息矩阵分别为

$ \begin{eqnarray} &&\frac{\partial l(\theta)}{\partial \theta} = \sum\limits_{i=1}^{N}(Y_{i}-p_{i})X_{i}, \notag \\ \label{eq5} &&-\frac{\partial^2 l(\theta)}{\partial \theta \partial \theta'} = \sum\limits_{i=1}^{N}p_{i}(1-p_{i})X_{i}X_{i}', \end{eqnarray} $

(2.5)

从而可得求解极大似然估计 $\hat\theta$的Newton-Raphson迭代公式为

$ \begin{eqnarray*}\label{eq6} \theta^{(n+1)} = \theta^{(n)} + \left(\sum\limits_{i=1}^{N}p_{i}^{(n)}(1-p_{i}^{(n)})X_{i}X_{i}'\right)^{-1} \left(\sum\limits_{i=1}^{N}(Y_{i}-p_{i}^{(n)})X_{i}\right), \end{eqnarray*} $

(2.6)

其中 $p_{i}^{(n)}=\frac{e^{X_i'\theta^{(n)}}}{1+e^{X_i'\theta^{(n)}}}$为 $p_i$的第 $n$次迭代值.

下面我们来应用QLB算法求解Logistic模型中回归参数的极大似然估计.我们知道, 算法关键是寻找二次下界矩阵 $B$, 使得 $B$是二阶导数矩阵 $\left(\frac{\partial^2 l(\theta)}{\partial \theta \partial \theta'}\right)$的一个下界. Böhning & Linday (1988)^[2]通过对(2.5) 式的观察得到: (1) $X_{i}X_{i}'$是非负定矩阵; (2) $p_{i}(1-p_{i})\le \frac{1}{4}, \ i=1, \cdots, N$.因此有

$\begin{eqnarray*} \frac{\partial^2 l(\theta)}{\partial \theta \partial \theta'} = -\sum\limits_{i=1}^{N}p_{i}(1-p_{i})X_{i}X_{i}' \geq -\frac{1}{4}\sum\limits_{i=1}^{N}X_{i}X_{i}', \end{eqnarray*} $

对任意 $\theta\in\Theta$均成立.从而找到了QLB算法中的一个二次下界矩阵 $B=-\frac{1}{4}\sum\limits_{i=1}^{N}X_{i}X_{i}'. $根据均值不等式可知, 这里的矩阵 $B$是二阶导数矩阵的最大下界, 也就是说此处的矩阵 $B$是二阶导数矩阵所有下界的上确界.因此, 我们可得求解 $\hat\theta$的QLB算法迭代公式如下:

$ \begin{eqnarray*}\label{eq7} \theta^{(n+1)} = \theta^{(n)} + \left(\frac{1}{4}\sum\limits_{i=1}^{N}X_{i}X_{i}'\right)^{-1} \left(\sum\limits_{i=1}^{N}(Y_{i}-p_{i}^{(n)})X_{i}\right). \end{eqnarray*} $

(2.7)

Böhning & Linday (1988) 证明了此迭代算法是单调线性收敛的^[2].由于QLB算法中的矩阵 $B$只与样本有关, 若矩阵 $\left(\sum\limits_{i=1}^{N}X_{i}X_{i}'\right)$不奇异或不接近奇异, 那么仅需要在整个迭代过程计算一次矩阵的逆即可.

进一步, 对于极大似然估计 $\hat\theta$方差的估计问题, 在Newton-Raphson算法下, 由于极大似然估计渐近正态性成立, 我们可以利用渐近方差的样本经验估计作为 $\hat\theta$方差的估计, 即 $\hat\theta$的协方差估计为 $\widehat{\text{Cov}}( \hat\theta ) = (-\frac{\partial^2 l(\hat{\theta})}{\partial \theta \partial \theta'})^{-1}$, 从而 $\hat\theta$的标准差(Standard Errors)为 $\widehat{\text{Cov}}( \hat\theta )$对角元素的平方根.而在QLB算法下, 我们则可以应用广泛使用的非参数Bootstrap重抽样方法^{[23, 24]}估计 $\hat\theta$的方差.

本节中, 我们将QLB算法应用于Logistic回归模型中参数的估计问题.首先, 我们进行数值模拟, 分别应用QLB算法与Newton-Raphson算法计算参数的极大似然估计.然后, 应用QLB算法分析了三个实际数据.所有分析均应用R软件.

我们考虑下述Logistic模型:

$ P(Y=1)=\frac{e^{\theta_0+\theta_1X_1+\theta_2X_2}}{1+e^{\theta_0+\theta_1X_1+\theta_2X_2}}. $

先从成功概率为0.2的Bernoulli分布生成协变量 $X_1$, 从标准正态分布生成 $X_2$, 从而二元响应变量 $Y$可由逆变换法产生.为了更好地研究模拟结果, 我们考虑不同的参数取值组合.分别取 $\theta_0=1$, $\theta_1=-0.5, \ 0, \ 0.5$, $\theta_2=-0.5, \ 0, \ 0.5$.分别取样本量 $N=300, \ 400$, 并设定QLB算法下Bootstrap次数为500次.

对于每种参数组合, 比较两种算法的参数估计结果:

(1) $\hat\theta_{N}$:应用Newton-Raphson算法, 基于迭代公式(2.6) 计算所得的极大似然估计;

(2) $\hat\theta_{Q}$:应用QLB算法, 基于迭代公式(2.7) 计算所得的极大似然估计.基于1000次模拟下, 我们得到估计的均值(Mean), 估计的样本标准差(SD), 标准差估计的均值(SE), 以及95%正态区间估计覆盖率(CP).模拟结果请见表 1. 表 1中结果表明, 在每种模拟设定下, $\theta_1$和 $\theta_2$的极大似然估计均为无偏估计. SE与SD充分接近, 表明提出的标准差的估计很好地估计了极大似然估计的样本标准差, 说明提出的标准差估计在实际应用中的合理性.区间估计覆盖率均接近95%, 表明在考虑的有限样本下, 极大似然估计的渐近正态性表现优良.进一步, 当样本量增大时, 估计的效率提高.数值模拟结果也表明, QLB算法与Newton-Raphson算法所得的结果一致, 这表明, QLB算法作为Newton-Raphson算法的合理替代, 能很好地应用于Logistic回归参数的估计问题.

实例1(急性淋巴细胞性白血病病人数据)

我们首先分析了一个50位急性淋巴细胞性白血病病人的数据, 数据来源于薛毅等^[25].在入院治疗时, 研究测得了病人外辕血中的细胞数(Cell), 单位为千个/mm $^{3}$; 淋巴结浸润等级(Lymph), 分为0, 1, 2, 3级; 出院后有无巩固治疗(Consolidation), 0表示无巩固治疗, 1表示有巩固治疗.因变量 $Y$为病人的生存时间, $Y=0$表示生存时间在1年以内, $Y=1$表示生存时间在1年或1年以上.

首先, 对淋巴结浸润等级(Lymph)以及出院后有无巩固治疗(Consolidation)变量的频数表(表 2)和条形图(图 1)分析发现, 出院后有无巩固治疗对病人生存时间影响显著, 即有巩固治疗的病人生存时间比没有巩固治疗的病人生存时间更长.而淋巴结浸润等级变量对病人的生存时间的影响没有明显的规律可循.

基于上述影响因素, 在如下Logistic模型框架下:

$ \log\frac{P(Y=1)}{1+P(Y=1)}=\theta_{0}+\theta_{1}\, {\hbox{Cell}}+\theta_{2}\, {\hbox{Lymph}}+\theta_{3}\, {\hbox{Consolidation}}, $

对数据进行了进一步的分析, 应用QLB算法计算回归参数的极大似然估计(Est), 估计的标准差(SE)仍由Bootstrap方法求得, Bootstrap次数设为1500次, 进而得到参数的95%置信区间(CI), 以及回归参数显著性检验的 $p$值( $p$-value).数据分析结果请见表 3.由表中结果可得, 出院后有无巩固治疗对病人生存时间有显著性影响(Est=2.8304, $p$-value＜0.001), 表明出院后接受巩固治疗能显著延长病人的生存时间.而外辕血中的细胞数与淋巴结浸润等级对病人生存时间无显著性影响.

实例2(教职工升职情况数据)

这里, 我们分析一个教职工升职情况的调查数据, 数据来源于王静龙等^[26].数据给出了某校5221位教职工升职情况以及性别(Sex)、工龄(Age)、文化程度(Education)三个相关影响因素.工龄取值为1, 2, 3, 4, 分别表示工龄小于5年, 6-15年, 16-29年以及大于等于30年; 文化程度取值为1, 2, 3, 4, 分别表示中学、大学、硕士、博士; 性别变量取值0, 1, 其中0表示女性, 1表示男性.因变量 $Y$取值0, 1, 其中0表示无升职, 1表示升职.通过频率表和条形图可以看出:随着工龄的增加, 或者教育程度的提高, 升职的人数逐渐增加(表 4, 图 2-3).因此可初步表明升职与否与工龄和教育程度有着明显相关性.

进一步, 建立升职与否与上述协变量之间的Logistic回归模型如下:

$ \log\frac{P(Y=1)}{1+P(Y=1)}=\theta_{0}+\theta_{1}\, {\hbox{Sex}}+\theta_{2}\, {\hbox{Age}}+\theta_{3}\, {\hbox{Education}}. $

应用QLB算法, 设定Bootstrap次数为500次, 所得数据分析结果请见表 5, 结果表明:工龄(Est=0.5053, $p$-value＜0.001) 与文化程度(Est=0.4385, $p$-value＜0.001) 对升职有显著性影响, 特别地, 工龄越高, 文化程度越高, 越有利于升职.结果同时也表明了性别对升职没有显著性影响.

实例3(农村居民健康行为数据)

最后, 我们讨论一个农村居民健康行为及其影响因素的研究数据, 数据来源于方积乾等^[27].在这项调查研究中, 随机调查了9760名15-55岁的居民, 其中吸烟男性有2229人, 过去1年内有过戒烟行为的有393人.此项研究的主要目的是分析和寻找对“戒烟”这一行为有影响的相关因素.研究中对多个指标进行合并得到了主要影响因素如下:年龄 $X_{1}$(岁), 家庭人均年实际收入 $X_{2}$(元), 婚姻状况 $X_{3}$(未婚=0, 已婚=1, 离婚丧偶=2), 职业 $X_{4}$(学生、家务及无业者=0, 农业劳动=1, 城乡农民工=2, 干部职工等其他=3), 受教育程度 $X_{5}$(文盲或半文盲=0, 小学=1, 初中=2, 高中及以上=3), 饮酒 $X_{6}$(不饮酒=0, 饮酒=1), 主动获取保健知识 $X_{7}$(不经常获取=0, 经常获取=1), 开始吸烟年龄 $X_{8}$(岁), 常去公共场所有无禁烟规定 $X_{9}$(无=0, 有=1), 家庭住室有无吸烟限制 $X_{10}$(无=0, 有=1), 有无经医生诊断的慢性病 $X_{11}$(无=0, 有=1), 自我健康总体评价 $X_{12}$(一般或较差=0, 较好=1, 很好=2).响应变量 $Y$表示过去1年内是否有过戒烟行为(未曾戒烟=0, 曾戒烟=1).进一步, 对原始数据中年龄 $X_{1}$以及开始吸烟年龄 $X_{8}$进行处理:对年龄 $X_{1}$, 年龄低于40岁的赋值为0, 大于40岁的赋值为1;对开始吸烟年龄 $X_{8}$, 小于25岁赋值为0, 大于25岁赋值为1.对收入 $X_{2}$, 做中心标准化处理.

研究上述所有潜在影响因素与戒烟与否之间的关系, 建立相应Logistic模型, 应用QLB算法并设定Bootstrap次数为1500次, 所得数据分析结果请见表 6.结果表明, 年龄超过40岁的男性更难戒烟, 而开始吸烟年龄大于25岁的男性更倾向于戒烟.受教育程度越高, 家庭住室有吸烟限制、经医生诊断患有慢性病, 以及自我健康总体评价较低的男性更倾向于戒烟.