离散选择模型中固定随机误差项方差的原因及后果

0. 引言

实证分析中，对于被解释变量为离散【注：如表示某结局是否发生的二值变量、表示数种选择方案的无序多分类变量、表示幸福程度的有序多分类变量（序数变量）等】的情况，研究者往往采用“离散选择模型（Discrete Choice Model）”对解释变量与被解释变量间的关系加以刻画。例如，最为我们所熟稔的Probit模型和Logistic模型便均为离散选择模型。然而，许多人或许未曾注意过，Probit模型和Logistic模型的参数识别需要满足一个隐含假设，即随机误差项的方差被固定为某一具体的值。基于此识别条件，随之而来地，是一些导致不利于我们分析和解读模型估计结果的问题。不失一般性，后文的讨论以被解释变量为二值的情形为例。

1. 线性概率模型

在涉及本文的主题之前，我们不妨先讨论一下对于二值被解释变量采用线性模型回归时的情况。考虑如下回归模型：

\(y_i=x_i\beta+\epsilon_i\)

因为\(y_i\)服从二点分布，其条件期望即被解释变量取值为1时的概率：

\(E(y_i|x_i)=P(y_i=1|x_i)=x_i\beta\)

此时回归系数反映的是解释变量的变化对于条件概率\(P(y_i=1|x_i)\)的边际影响，且该影响是恒定的，因此上述模型也被称作“线性概率模型（Linear Probability Model）”。例如，假设\(y_i\)表示是否购买养老保险，取1时为购买而取0时为不购买，第\(j\)个解释变量\(x_{ij}\)为年龄，边际影响\(\frac{\partial{P(y_i=1|x_i)}}{\partial{x_{ij}}}=0.038\)意味着控制其他变量不变的条件下，年龄增大10岁，购买养老保险的概率提高0.38。尽管线性概率模型相对易于解释，且非常容易通过OLS加以估计，但它存在着三点不足。

1.1 随机误差项的非正态性

显而易见，由于被解释变量服从二点分布，随机误差项\(\epsilon_i\)同样服从二点分布，以\(1-F(x_i\beta)\)和\(F(x_i\beta)\)概率取值为\(-x_i\beta\)和\(1-x_i\beta\)。不过，\(\epsilon_i\)违反正态性假设不会对线性概率模型的估计产生影响，因为有限样本下OLS估计量的无偏性并不要求\(\epsilon_i\)服从正态分布，而在大样本的情况下，OLS估计量的分布将收敛于正态分布，从而也不会对统计推断造成困扰。

1.2 随机误差项的异方差性

基于与上述相同的理由，\(\epsilon_i\)是（条件）异方差的，即\(Var(\epsilon_i|x_i)=x_i\beta(1-x_i\beta)\)依赖于\(x_i\)。尽管违反（条件）同方差假设不会影响OLS估计量的无偏性、一致性和渐近正态性，但OLS估计量将不再是（渐近）有效的。不过，该问题可以通过加权最小二乘或者采用Eicker-Huber-White标准误等方法解决。

1.3 不满足\(0\leq E(y_i|x_i)\leq 1\)

随机误差项的非正态性和异方差性对于估计线性概率模型来说其实属于不是问题的问题，因为我们总有方法加以克服。线性概率模型最大的缺点在于，我们知道概率必须介于0和1之间，而我们无法将\(x_i\beta\)限定在这一区间内，从而导致出现毫无意义的概率和负的方差。

2. Probit模型与Logistic模型的识别

既然应用线性概率模型的最大问题在于不满足\(0\leq E(y_i|x_i)\leq 1\)这一约束条件，那末一个自然而然的想法是通过一种技术来保证\(y_i\)的条件概率一定位于0和1之间。注意1.1小节中的细节，\(F(x_i\beta)=P(y_i=1|x_i)\)，\(F(·)\)为逆连接函数（Inverse Link Function）【注：此处的表述在机器学习的语境下略有出入，可参见笔者的另一篇文章《[CuPy for Machine Learning] Logistic分类器》】。在线性概率模型中，\(F(·)\)为恒等（Identity）连接，所以\(P(y_i=1|x_i)\)的值域为全体实数。那末，如果逆连接函数为分布函数，则可以保证\(0\leq E(y_i|x_i)\leq 1\)。为此，考虑标准正态分布的分布函数作为逆连接函数：

\(P(y_i=1|x_i)=F(x_i\beta)=\Phi(x_i\beta)=\int_{-\infty}^{x_i\beta}\phi(z)dz\)

即得到Probit模型。

或者选择Logistic分布的分布函数：

\(P(y_i=1|x_i)=F(x_i\beta)=\Lambda(x_i\beta)=\frac{e^{x_i\beta}}{1+e^{x_i\beta}}\)

得到Logistic模型。

然而，若要说明Probit模型和Logistic模型的识别问题，即本文开篇所说的“随机误差项的方差被固定为某一具体的值”，仅仅依靠上述推导是不够的。我们必须借助潜变量模型（Latent Variable Model）：

\(y_i^*=x_i\beta+\epsilon_i\)

式中，\(y_i^*\)为不可观测的潜变量，随机误差项服从某一关于0的对称分布\(F\)，即\(\epsilon_i\sim F(\epsilon_i)\)。我们实际观测到的\(y_i\)为：

\(y_i=\begin{cases}1,y_i^*>0\\0,y_i^*\leq 0\end{cases}\)

因此

\(\begin{aligned}E(y_i|x_i)&=P(y_i=1|x_i)\\&=P(y_i^*>0|x_i)\\&=P(\epsilon_i>-x_i\beta|x_i)\\&=P(\epsilon_i<x_i\beta|x_i)\\&=F(x_i\beta)\end{aligned}\)

如果\(F\)为标准正态分布，则得到Probit模型；如果\(F\)为Logistic分布，则得到Logistic模型。

现在，我们假设\(e_i=\sigma \epsilon_i\)，其中\(\epsilon_i\)服从一个具有单位方差的分布\(F\)，则：

\(\begin{aligned}E(y_i|x_i)&=P(y_i=1|x_i)\\&=P(y_i^*>0|x_i)\\&=P(e_i>-x_i\beta|x_i)\\&=P(e_i<x_i\beta|x_i)\\&=F(x_i\frac{\beta}{\sigma})\\&=F(x_i\beta^*)\end{aligned}\)

显而易见，此时我们无法同时识别\(\beta\)和\(\sigma\)而只能识别二者之比\(\beta^*=\frac{\beta}{\sigma}\)【注：模型的识别是一个设定问题而不是估计问题。例如，结构方程模型中要求内生潜变量间的结构系数矩阵\(B\)满足\((I-B)\)为非奇异矩阵，否则模型无法识别；又例如，需要将潜变量的方差固定为1或者将某一显变量的因子载荷固定为常数（一般为1），否则因子载荷和潜变量的方差便无法识别】。因此，为了识别系数\(\beta\)，必须将\(\sigma\)设定为固定的值，通常标准正态分布为1而Logistic分布为\(\pi/\sqrt{3}\)。在一定意义上，可以将\(\sigma\)理解为重新度量了模型系数\(\beta\)的尺度因子（Scaling Factor）。当然，固定\(\sigma\)的值并不会影响到我们对\(P(y_i=1|x_i)\)的解释，而且事实上完全可以为\(\sigma\)设定任意的值，将之设定为1或者\(\pi/\sqrt{3}\)仅仅是因为它们与单位方差相同【注：Long(1997)对此作了一个简要的证明，同时他也阐述了为何我们偏好于选择标准正态分布和Logistic分布的分布函数作为连接函数】。

3. 模型的系数比较问题

尽管固定\(\sigma\)是Probit模型或者Logistic模型可识别的必要条件，但是这也导致了一个实际操作中不可避免的问题，即无法比较不同模型之间系数的大小。原因在于随着纳入模型的变量增加或减少，\(\sigma\)也随之改变，而作为尺度因子，\(\sigma\)的改变必然引起系数估计值的改变。出于模型识别的要求，\(\sigma\)被设为固定值，因此无法捕捉模型的变异。为直观地说明该问题，考虑如下潜变量模型【注：此处以Probit模型为例，但结论也同样适用于Logistic模型】：

\(y_i^*=\beta_0+\beta_1x_{1i}+\beta_2x_{2i}+\sigma\epsilon_i\)

其中\(\sigma=1\)，\(\epsilon_i\sim N(0,1)\)。如果\(x_{1i}\)与\(x_{2i}\)相关，但是估计时\(x_{2i}\)被遗漏，则导致遗漏变量偏差。由于遗漏变量偏差的情况也同样出现在线性回归模型中，所以不再赘述。我们重点关注\(x_{1i}\)与\(x_{2i}\)不相关的情况。估计上式，有：

\(\frac{\beta_1}{\sigma}=\beta_1^F\)

现在从模型中剔除\(x_{2i}\)，则误差方差变为：

\(\beta_2^2Var(x_{i2})+Var(\epsilon_i)\)

此时\(\sigma\)变为：

\(\sigma=\sqrt{\beta_2^2Var(x_{2i})+1}\)

\(x_{1i}\)的系数相应发生改变：

\(\frac{\beta_1}{\sqrt{\beta_2^2Var(x_{2i})+1}}=\beta_1^S\neq \beta_1^F\)

因为\(\sqrt{\beta_2^2Var(x_{2i})+1}>1\)，除非\(\beta_2=0\)或者\(Var(x_{2i})=0\)，否则总有\(|\beta_1^S|<|\beta_1^F|\)。如果我们将同时估计\(x_{1i}\)和\(x_{2i}\)的模型视作全模型（用上标\(F\)表示），仅估计\(x_{1i}\)的模型视作子模型（用上标\(S\)表示）。因为存在着“衰减偏差（Attenuation Bias）”，显而易见，我们无法直接比较嵌套模型中同一变量系数的大小。因此，针对Probit（Logistic）模型系数不可比的问题，可以采用Karlson-Holm-Breen分解（Karlson et al.，2013）、计算平均偏效应（Wooldridge，2010）或者计算\(y_i^*\)标准化系数（Long，1997）等方法加以解决。当然，如果样本拟合值中小于0或大于1的数据不多，那末估计前文所介绍的线性概率模型也未尝不可。

4. 模拟

通过Stata模拟来验证上述讨论。以Probit模型为例，结果如下图所示。可以看到，子模型\(x_{i1}\)的系数较之全模型明显减小，而两模型在样本均值处计算的边际效应也存在较大差别【注：关于这一问题的讨论，参见Mood（2010）】。然而，平均偏效应和\(y_i^*\)标准化系数则十分接近，意味着不同模型系数比较时可以选择汇报平均偏效应或者\(y_i^*\)标准化系数。

相关Stata代码如下：

drop _all
set obs 10000
set seed 196
tempname add

/* Generate random data. */

gen e = rnormal()
gen x1 = rnormal()
gen x2 = rnormal()
gen ystar = x1 + 2*x2 + e
gen y = (ystar > 0)

/* Estimation of full model. */

qui probit y x1 x2, r
qui margin, dydx(x1)
mat `add' = r(b)
qui estadd scalar ape = `add'[1,1]
qui margin, dydx(x1) atmeans
mat `add' = r(b)
qui estadd scalar mfx = `add'[1,1]
qui listcoef, std constantoff
mat `add' = r(table)
qui estadd scalar ystd = `add'[1,5]
eststo m1

/* Estimation of reduced model. */

qui probit y x1, r
qui margin, dydx(x1)
mat `add' = r(b)
qui estadd scalar ape = `add'[1,1]
qui margin, dydx(x1) atmeans
mat `add' = r(b)
qui estadd scalar mfx = `add'[1,1]
qui listcoef, std constantoff
mat `add' = r(table)
qui estadd scalar ystd = `add'[1,5]
eststo m2

/* Display the results. */

esttab m1 m2, not compress ///
	mtitle("Model A" "Model B") ///
	title("Comparison of estimates (N=10000)") ///
	nonumbers nodepvars drop(x2 _cons) ///
	stats(ape mfx ystd) ///
	addnotes("Model A with control for x2.")

5. 广义线性模型

最后，我们来讨论一下如何在更一般化的框架下理解Probit模型或Logistic模型，即广义线性模型（Generalized Linear Model）的框架。广义线性模型是一般线性模型的推广，允许基于被解释变量的性质更加灵活地设定模型，同时又通过转化线性模型的系统成分而维持了一般线性模型的结构。广义线性模型主要由两个部分组成，误差结构和连接函数。前者要求被解释变量的方差为其均值的函数，因此被解释变量服从某一特定的指数族分布是应用广义线性模型的前提；后者是被解释变量均值的转化，转化的结果为回归参数的线性函数，在广义线性模型的语境中被称作线性预测（Linear Predictor），代表模型的系统成分。

任一随机变量\(Y\)的概率密度函数属于指数族的前提是它可以被转化为如下的典型形式【注：不同文献中关于指数族分布的函数形式可能有所出入，但是本质相同】：

\(f(y|\theta)=e^{y\theta+b(\theta)+c(y)}\)

上式中\(\theta\)被称作自然参数（Natural Parameter），它是分布均值的一个函数——即\(\theta(\mu)\)。此时连接函数\(g(\mu)=\theta=\eta=X\beta\)，因此\(\theta\)也是任意指数族分布的“标准连接函数（Canonical Link Function）”。\(b(\theta)\)被称作累积量函数（Cumulant Function），用于描述分布的矩（Moment）。\(c(y)\)被称作正则项（Normalization Term），用于保证概率密度的正则性【注：\(\int_{-\infty}^{\infty}{f(x)}dx=1\)】。

回到本文的主题上，由于被解释变量\(y_i\)服从二点分布，即\(y_i\sim B(\pi)\)。二点分布的概率密度（质量）函数为：

\(f(y_i|\pi)=\pi^{y_i}(1-\pi)^{1-y_i}\)

将之转化为指数族分布：

\(\begin{aligned}f(y_i|\pi)&=\pi^{y_i}(1-\pi)^{1-y_i}\\&=e^{y_i\ln\pi+(1-y_i)\ln{1-\pi}}\\&=e^{y_i\ln{\frac{\pi}{1-\pi}}+\ln{1-\pi}}\\&=e^{y_i\theta-\ln{1+e^\theta}}\end{aligned}\)

其中自然参数\(\theta=\ln{\frac{\pi}{1-\pi}}\)即标准连接函数。因为\(\pi=E(y_i|x_i)=P(y_i=1|x_i)\)，所以\(\ln{\frac{P(y_i=1|x_i)}{1-P(y_i=1|x_i)}}=x_i\beta\)。进一步地，\(P(y_i=1|x_i)=g^{-1}(\theta)=\frac{e^{x_i\beta}}{1-e^{x_i\beta}}\)。如此便得到Logistic模型。

事实上，关于连接函数的选择是相对灵活的。除标准连接函数外，也可以选择非标准连接函数（Noncanonical Link Function）【注：只要该函数是单调可微的】。尽管标准连接函数具有良好的统计性质【注：如\(x'_iy_i\)是关于参数\(\beta\)的充分统计量，或者线性预测\(\eta\)的结果被严格映射在分布均值的取值范围内等。此外，从参数估计的角度来看，采用标准连接函数时对数似然函数的观测海森矩阵等价于其期望海森矩阵，因此Newton-Raphson法等价于迭代加权最小二乘法（Iteratively Reweighted Least Squares，IRLS）】，但是在某些情形下，非标准连接函数的模型拟合效果可能优于标准连接函数，同时也令关于模型系数的解释更加灵活。例如，针对二点分布可以选择\(g(\pi)=ln(\pi)\)作为连接函数，以便将系数解释为相对风险（Relative Risks）而不是几率比（Odds Ratio），此时模型被称作Log-binomial模型，常见于流行病学研究中；而Probit模型则是以\(g(\pi)=\Phi^{-1}(\pi)\)作为连接函数。

参考文献

[1]Cameron A C, Trivedi P K. Microeconometrics: Methods and applications[M]. Cambridge University Press, 2005.

[2]Hardin J W, Hilbe J M. Generalized estimating equations[M]. CRC Press, 2003.

[3]Karlson K B, Holm A, Breen R. Comparing regression coefficients between same-sample nested models using Logit and Probit: A new method[J]. Sociological Methodology, 2013, 42(1): 286-313.

[4]Mood C. Logistic regression: Why we cannot do what we think we can do, and what we can do about it[J]. European Sociological Review, 2010, 26(1): 67-82.

[5]Long J S. Regression models for categorical and limited dependent variables[M]. SAGE Publications, 1997.

[6]Wooldridge J M. Econometric analysis of cross selection and panel data (2e.)[M]. MIT Press, 2010.

[创作不易，转载请注明出处]