Firth惩罚逻辑回归 | Codeboy's Blog

当二元数据被一条超平面分开时，原似然会鼓励系数不断增大。Firth方法给似然增加一项随信息矩阵变化的惩罚，使有限数据也能产生有限估计。要使用这一结果，应把三个问题分开：目标为什么有有限解，求解器怎样逼近它，所报区间是否具有要求的覆盖率。

形式陈述

模型和被优化的目标

给定满列秩 $X\in\mathbb R^{n\times d}$ ，n≥d≥1，条件于设计的响应相互独立，且

Y_i\sim\operatorname{Bernoulli}(p_i),\qquad p_i=\frac{e^{x_i^\mathsf T\beta}}{1+e^{x_i^\mathsf T\beta}},\qquad \beta\in\mathbb R^d.

这是广义线性模型的Bernoulli–logit情形。记观测对数似然与Fisher信息为

\ell(\beta)=\sum_i\{y_i x_i^\mathsf T\beta-\log(1+e^{x_i^\mathsf T\beta})\},\qquad I(\beta)=X^\mathsf T W X,\qquad w_i=p_i(1-p_i).

Firth惩罚逻辑回归选择下面目标的全局极大点：

\widetilde\ell(\beta)=\ell(\beta)+\tfrac12\log\det I(\beta).\tag{1}

有限β下I正定，其行列式严格正，因此式(1)有定义。在这个满列秩模型中，式(1)至少有一个有限全局极大点，无论原资料是否分离。此结论不附带全局唯一性或任意迭代法的收敛保证。

这里限定logit链接和β坐标。一般模型也可以加入Jeffreys信息行列式惩罚，但它未必等于Firth的一阶均值偏差修正；更换非线性参数坐标后，密度的极大点也未必保持不变。

修正得分的逐项推导

令

h_i=w_i x_i^\mathsf T I^{-1}x_i.

它们是矩阵 $W^{1/2}X(X^\mathsf T W X)^{-1}X^\mathsf T W^{1/2}$ 的对角元。该矩阵是秩d的正交投影，所以 $0\le h_i\le1$ 、 $\sum_i h_i=d$ 。

对可逆矩阵，沿方向的行列式微分为 $d\log\det I=\operatorname{tr}(I^{-1}dI)$ ；可由 $\det(I+tA)=1+t\operatorname{tr}A+O(t^2)$ 得到。又有

\nabla_\beta w_i=w_i(1-2p_i)x_i.

把 $dI=\sum_i(dw_i)x_ix_i^\mathsf T$ 代入并取迹，得到

\nabla\widetilde\ell(\beta) =U^*(\beta) =X^\mathsf T\{y-p+h\mathbin{\odot}(\tfrac12\mathbf1-p)\}.\tag{2}

符号⊙表示逐分量相乘。修正量依赖当前β及整个设计，不能把一般回归解释成预先给每条记录固定增加半个成功。任一有限极大点必须使式(2)为零；反过来，零得分只是驻点条件。

为什么参数不会逃向无穷

将行列式按多线性展开，有有限子式恒等式

\det(X^\mathsf T W X) =\sum_{S\subseteq\{1,\ldots,n\},\ |S|=d} \det(X_S)^2\prod_{i\in S}w_i.\tag{3}

这里X_S由选中的d行组成。为看出式(3)，令A=W^{1/2}X，展开 $\det(A^\mathsf T A)$ ：重复选同一行的项因交替性消失；每组不同的d行合起来给 $\det(A_S)^2$ 。这就是所需的Cauchy–Binet恒等式。

只需考虑 $\det X_S\ne0$ 的子集。X_S可逆，所以存在 $c_S>0$ ，使

\max_{i\in S}|x_i^\mathsf T\beta|\ge c_S\|\beta\|_2

对所有β成立。例如可取 $c_S=1/(\sqrt d\|X_S^{-1}\|_2)$ 。另一方面，写η=x_i^Tβ，有

0<w_i=\frac{e^{-|\eta|}}{(1+e^{-|\eta|})^2}\le e^{-|\eta|},\qquad w_i\le\tfrac14.

故式(3)中每个非零子式项至多是一个常数乘 $e^{-c_S\|\beta\|_2}$ 。非零子式有限且至少有一个，取它们的最小c_S，便有常数A,c>0满足

\det I(\beta)\le A e^{-c\|\beta\|_2}.

原Bernoulli对数似然始终≤0，因此 $\widetilde\ell(\beta)\le\tfrac12\log A-\tfrac c2\|\beta\|_2\to-\infty$ 。足够大球外的值比 $\widetilde\ell(0)$ 小；在闭球内用极值定理，就取得有限全局极大点。证明控制所有逃逸方向，而非仅检查固定的一条直线。

直觉

靠把概率推到0或1来提高拟合时，部分信息权重 $p_i(1-p_i)$ 会变得很小。惩罚项监测整个信息椭球的体积：只要参数逃得足够远，每个能确定全部系数的满秩子设计都会失去一些信息，行列式最终趋零。

修正得分中的 $h_i(1/2-p_i)$ 则把这种整体变化写成每行的贡献。高杠杆行影响较大，当前概率位于1/2哪一侧决定修正的符号。它不是固定的欧氏长度惩罚，也不保证每一个回归系数的绝对值都缩小。

“偏差减少”指指定坐标和渐近模型下的偏差阶数，不等于每一份数据上更接近真值。下面的单截距计算既给出阶数抵消，也能直接看到有限样本下仍然有偏。

例子与边界

单截距：把端点变为有限值

若X只有一列1，令 $S=\sum_iY_i$ ，p=σ(β)。忽略不依赖β的常数，式(1)为

\widetilde\ell(\beta)=(S+\tfrac12)\log p+(n-S+\tfrac12)\log(1-p).

对β求导得 $S+1/2-(n+1)p$ ，故唯一极大点是

\widetilde p=\frac{S+1/2}{n+1},\qquad \widetilde\beta=\log\frac{S+1/2}{n-S+1/2}.\tag{4}

S=0或n时仍有限。例如n=1，两种可能估计为−log3与log3。若真p=3/4，则 $E\widetilde\beta=\tfrac12\log3$ ，真β=log3，偏差并没有变成零。

单截距的偏差阶数：抵消发生在哪里

固定真值 $p\in(0,1)$ ，令q=1−p，n趋无穷， $S\sim\operatorname{Binomial}(n,p)$ 。设 $a=1/2-p$ 、 $\Delta=\widetilde p-p=(S-np+a)/(n+1)$ ，则

E\Delta=\frac a{n+1},\qquad E\Delta^2=\frac{npq+a^2}{(n+1)^2}.

对 $f(u)=\log\{u/(1-u)\}$ ，有 $f'(p)=1/(pq)$ 、 $f''(p)=(2p-1)/(p^2q^2)$ 。应用Taylor展开时，一次与二次项的期望之和满足

f'(p)E\Delta+\tfrac12f''(p)E\Delta^2 =\frac1n\left\{\frac{1/2-p}{pq}+\frac{2p-1}{2pq}\right\}+O(n^{-2}) =O(n^{-2}).\tag{5}

还需控制余项，不能只看到括号为零便结束。令V=S−np，由独立Bernoulli中心矩展开，

EV^2=npq,\qquad EV^3=npq(1-2p),\qquad EV^4=3(npq)^2+npq(1-6pq).

三阶矩中只留下同一索引，四阶矩中只留下四次同一索引或两对索引，因而可直接得到这些式子。代入Δ后， $E\Delta^3=O(n^{-2})$ 、 $E\Delta^4=O(n^{-2})$ 。

取固定ε>0使 $[p-\varepsilon,p+\varepsilon]\subset(0,1)$ 。在 $|\Delta|\le\varepsilon$ 上，展开至三次，四阶导数有界，故余项绝对期望为O(n^{-2})。在补集上，Hoeffding界给 $P(|\Delta|>\varepsilon)\le2e^{-c_p n}$ ，对所有足够大n成立；而式(4)保证 $|\widetilde\beta|\le\log(2n+1)$ ，Δ也有界。因此尾部对展开及真实函数的贡献至多为 $O((1+\log n)e^{-c_p n})$ 。综合式(5)，

E_p\widetilde\beta-\log\frac p{1-p}=O(n^{-2}).\tag{6}

这是固定内部p的结论，常数可以依赖p，不是对趋向0或1的参数序列给统一保证。普通无约束logit MLE在S=0、n时没有有限值，这两种事件每个有限n都有正概率；不能把其未定义的无条件均值直接写成有限的偏差来相减。

Jeffreys测度与密度极大点的坐标区别

在β坐标下，Jeffreys因子正比 $\sqrt{np(1-p)}$ ，乘似然后的β密度极大点给式(4)。变到p坐标时，要连同Jacobian $d\beta/dp=1/[p(1-p)]$ 一起改变密度，后验密度便正比

p^{S-1/2}(1-p)^{n-S-1/2}.

它是Beta(S+1/2,n−S+1/2)密度；在S>1/2且n−S>1/2时，其内部众数是 $(S-1/2)/(n-1)$ ，通常不同于式(4)。Jeffreys测度的不变性没有使非线性坐标下的众数也不变。这个Beta后验的均值恰好等于式(4)，只是本例的代数一致，不能据此把一般Firth回归叫成后验均值法。

可逆线性换坐标β=Aγ则不同：信息行列式只多出常数因子 $\det(A)^2$ ，惩罚目标只差常数，拟合概率保持不变。

满秩是有限性证明的必要输入

若设计秩亏，I的行列式恒为0，式(1)不是一个可用于拟合的有限实目标。先删除不可识别方向或加入明确的识别约束，再在满秩坐标下计算。把零行列式替换成某个软件小常数，会另行改变目标，并不自动继承本页定理。

有限系数也没有承诺有限样本的Wald覆盖。n=1时，式(4)给 $\widetilde p=1/4$ 或3/4，所以 $I(\widetilde\beta)=3/16$ 。用该信息计算标准误为 $4/\sqrt3$ ；两种名义95%Wald区间的上端都不超过 $\log3+1.96(4/\sqrt3)<5.626$ 。若真β=6，两种区间都会漏掉它，实际覆盖率为0。这不违反有限性，也不违反固定内部p、n→∞的式(6)。

推论与应用

一次可检查的求解更新

从有限β开始，计算p、w、I、h与U*，解线性方程

I\delta=U^*.

若U*≠0，则 $(U^*)^\mathsf T\delta=(U^*)^\mathsf TI^{-1}U^*>0$ ，δ是惩罚目标的上升方向。选择c∈(0,1)、ρ∈(0,1)，从s=1开始缩小s←ρs，直到

\widetilde\ell(\beta+s\delta)\ge\widetilde\ell(\beta)+cs(U^*)^\mathsf T\delta.

目标在有限点光滑，正方向导数保证足够小s会通过。更新后重新计算h，不能将上一步的修正响应固定使用到底。

稠密设计的一轮信息构造、因子分解及所有h_i计算需O(nd²+d³)次算术；通过分解求解各 $I^{-1}x_i$ ，无需显式存储完整逆矩阵。每次回溯重新评价目标还需要相应计算。用稳定的log-sum-exp计算原似然，并从正定分解计算logdet，避免直接将极小权重相乘。

小梯度和小步长是局部停机诊断。报告还应保留目标值、信息矩阵条件与初始值敏感性；本页没有以一次停机标志证明全局最优。分离检查仍可用来解释为什么原MLE不存在，不能因为惩罚拟合成功就省略对原问题的说明。

参考资料

David Firth，Bias Reduction of Maximum Likelihood Estimates，Biometrika 80(1),1993,pp.27–38；偏差修正方法的原始论文。
Ioannis Kosmidis and David Firth，Jeffreys-prior Penalty, Finiteness and Shrinkage in Binomial-response Generalized Linear Models，Biometrika 108(1),2021,pp.71–82，§2.2、§4与补充S2；有限性、修正得分与推断边界。本文采用子式展开给出有限性的完整证明，并在单截距模型中直接证明式(6)。

形式陈述 ​

模型和被优化的目标 ​

修正得分的逐项推导 ​

为什么参数不会逃向无穷 ​

单截距：把端点变为有限值 ​

单截距的偏差阶数：抵消发生在哪里 ​

Jeffreys测度与密度极大点的坐标区别 ​

满秩是有限性证明的必要输入 ​

一次可检查的求解更新 ​

反向链接