feologit：固定效应有序Logit模型

发布时间：2021-01-07 阅读 19930

Stata连享会主页 || 视频 || 推文 || 知乎

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者： 庄子安 (中山大学)
邮箱： 1484712416@qq.com

指导老师： 连玉君 (中山大学，arlionn@163.com)

1. 应用背景

1.1 模型使用场景

你可能听说过固定效应模型，你也可能听说过有序 Logit 模型，那你有听说过固定效应有序 Logit 模型吗？本文将详细介绍固定效应有序 Logit 模型，以及该模型在 Stata 中的命令 feologit 的使用方法。

在多元 Logit 模型中，如果各个类别变量是有序的，但各个类别变量之间的差距未知，此时应该使用 有序 Logit 模型。

例如，在做性格测试问卷调查时，选项通常为「非常同意 (SA)」, 「同意 (A)」, 「不同意 (D)」和「非常不同意 (SD)」。「非常同意 (SA)」到「同意 (A)」的距离不一定等于「同意 (A)」到「不同意 (D)」的距离：

y_{i} = {\begin{array}{rcl} S D & τ_{0} = - \infty \leq y_{i}^{*} \leq τ_{1} \\ D & τ_{1} \leq y_{i}^{*} \leq τ_{2} \\ A & τ_{2} \leq y_{i}^{*} \leq τ_{3} \\ S A & τ_{3} \leq y_{i}^{*} \leq τ_{4} = + \infty \end{array}

同理，教育程度也是一个多类别变量，包括「小学」, 「初中」, 「高中」, 「本科及以上」。同样，我们无法界定不同学历之间的差距是否相同。

另外，如果数据为 面板数据，为了研究每个个体难以观测的不随时间变化的差异，例如个人消费习惯、企业文化、国家社会制度等，应结合采用固定效应模型。

因此，当数据特征为面板数据且被解释变量为有序类别变量时，应该采用 固定效应有序 Logit 模型 ：

panel data + ordered response = Fiexed effect ordered logit model

1.2 典型应用领域

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。

A. Anna Cristina D’Addio (2007): 工作满意度 (Job Satisfication) 的影响因素研究

该文献研究了当个体之间对工作满意度水平度量可能存在差异时，影响工作满意度的因素，即在前人对工作满意度影响因子的基础上加入了个体效应。

作者使用的数据源于 1995-1999 年欧洲共同体住户小组(ECHP)调查。这些数据按年收集，涉及家庭结构、家庭和家庭成员的收入和就业等几个问题，并提供了关于社会变化和个人行为动态的独特信息。

被采访者需要对自己的工作打一个分，范围为 1-6，作为实证模型的被解释变量。被解释变量有序但类别间的「距离」无法界定，且为面板数据，因此可以使用 固定效应有序 Logit 模型。

类似的文献还有：Anna Cristina D’Addio (2007)，Liliana Winkelmann and Rainer Winkelmann (1998)，Ada Ferrer-i-Carbonell and Paul Frijters(2004)。

B. Anne Boring (2014)：学生对教师评价中的性别偏见 (Gender Biases)

该文旨在研究 SET(Student Evaluation of Teachers) 中的性格差异问题，即不同性别的学生在给不同性别的老师做评价时是否会受到性别的外生影响。

作者使用了一个来自法国大学的独特数据库来研究学生对教师的评价存在性别偏见。从 2008 年起，学生需每年完成 SET 测评，该测评要求学生从四个维度评价教师：课程内容、作业和测验、授课风格、课程广度。学生对此需分别打出一个范围为0-4的分数，分数越高代表越满意。实证中，为了研究性别差异，最终使用了 固定效应有序 Logit 模型 。

2. 固定效应有序 Logit 模型

2.1 模型设定

固定效应有序 logit 模型使用潜在变量 $y^{*}$ 将可观测特征 $x$ 与可观测有序因变量 $y$ 相关联，可观测有序因变量 $y$ 可以取值 0-K。个体 $i$ 在时间 $t$ 的潜在变量 $y_{i t}^{*}$ 线性决定于 $x_{i t}$ 和两个无法观测的变量 $α_{i}$ 和 $ε_{i t}$ ：

y_{i t}^{*} = x_{i t}^{^{'}} β + α_{i} + ε_{i t} i = 1, \dots, N t = 1, \dots, T

其中， $α_{i}$ 为个体异质性截距项，且统计上取决于 $x_{i t}$ ， $ε_{i t}$ 为残差项。另外，潜在变量 $y_{i t}^{*}$ 与 $y_{i t}$ 的联系定义如下：

y_{i t} = k i f τ_{i k} < y_{i t}^{*} < τ_{i k + 1} k = 1, \dots, K

其中， $τ_{i k}$ 为不同个体 $i$ 的阈值。在固定效应有序 Logit 模型中，阈值可以因人而异。除了规定最低和最高阈值为负无穷大和正无穷大之外，关于个体特定阈值的唯一假设为，每个个体的阈值一直在增加：

τ_{i 1} = - \infty; - \infty < τ_{i k} < τ_{i k + 1} < \infty, \forall k = 2, \dots, K - 1; τ_{i K + 1} = \infty (1)

另外，该模型还假设残差项 $ε_{i t}$ 独立同分布于标准 Logistic 分布，因此 $ε_{i t}$ 的分布函数为：

F (ε_{i t} | x_{i t}, α_{i}) = F (ε_{i t}) = \frac{1}{1 + \exp (- ε_{i t})} \equiv Λ (ε_{i t})

因此，个体 $i$ 在时间 $t$ 的观测值等于 $k$ 的概率为：

Pr (y_{i t} = k | x_{i t}, α_{i}) = Λ (τ_{i k + 1} - x_{i t}^{'} β - α_{i}) - Λ (τ_{i k} - x_{i t}^{'} β - α_{i}) (2)

由等式 (2) 可知，概率不仅取决于 $x_{i t}$ 和 $β$ ，还取决于 $α_{i}$ ， $τ_{i}$ 和 $τ_{i k + 1}$ 。在有序 Logit 模型中，采用极大似然估计来得到 $β$ 和 $τ_{i k}$ ，但在该模型中，由于个体异质性截距项 $α_{i}$ 的存在，我们只能识别出 $τ_{i k} - α_{i} = α_{i k}$ ，产生了 Incidental Parameter Problem (伴随参数问题，Chamberlain 1980)，导致无法得到 $β$ 的一致估计量。解决方法以条件极大似然估 CML 为基础。

2.2 估计方法

2.2.1 CML estimator

固定效应有序 Logit 模型使用的估计方法是以 CML estimator 为基础的。在 Stata 中，clogit 就是使用该估计方法。类似的，feologit 也是以 clogit 为基础的。其原理很简单，就是通过条件极大似然估计摆脱 $α_{i}$ 的影响。具体估计过程如下：

定义二元变量 $d_{i t}^{k}$ 为在截断点 $k$ 处的有序变量二分所得： $d_{i t}^{k} = 1 (y_{i t} \geq k)$ 。我们假设个体 $i$ 的因变量被观察了 $T$ 次，定义：

g_{i}^{k} = \sum_{t = 1}^{T} d_{i t}^{k}

那么在 $g_{i}^{k} = \sum_{t = 1}^{T} d_{i t}^{k}$ 的条件下，观察到 $d_{i}^{k} = (d_{i 1}^{k}, \dots, d_{i T}^{k})^{'}$ 的概率为：

$P_{i}^{k} (β) = Pr (d_{i}^{k} | \sum_{t = 1}^{T} d_{i t}^{k} = g_{i}^{k}) = \frac{\exp (d_{i}^{k^{'}} x_{i} β)}{\sum_{j \in B_{i}} \exp ({j^{'} x}_{i} β)} (3)$

其中， $j$ 代表一个由 $T$ 个等于 $0$ 或 $1$ 的 $j_{t}$ 组成的 $T$ 维向量，且满足 $\sum_{t = 1}^{T} j_{i t}^{k} = g_{i}^{k}$ 。另外， $B_{i}$ 代表所有满足条件的集合，里面有 $g_{i}^{k}$ 个 $1$ 和 $T - g_{i}^{k}$ 个 $0$ 。

可见， (3) 式中的条件概率不取决于 $α_{i}$ 和阈值。因此可以得到条件对数似然函数：

{L L}^{k} (b) = \sum_{i = 1}^{N} \log P_{i}^{k} (b) (4)

通过极大化 (4) 式即可得到 $β$ 的一致估计量。但这种方法只使用了一个截断点，很多时候我们需要同时考察多个截断点。

2.2.2 BUC estimator

BUC estimator (Baetschmann, Staub, and Winkelmann [2015])结合了使用不同截断点 $k$ 的对数似然函数，其方程如下：

L L^{B U C} (b) = \sum_{k = 2}^{K} L L^{k} (b) (5)

其中， $L L^{k} (b)$ 来源于公式 $(4)$ 。BUC estimator 原理为最大化 $(5)$ 式，我们也可以把它理解为一个施加了约束条件的 CML estimator，因为它隐含了如下的约束条件： ${\hat{β}}^{2} = . . . = {\hat{β}}^{K}$ 。其原理分为三步：

样本中每个个体的观察值都被替换为 $K - 1$ 个自身的副本
每个副本在不同的截断点二分
使用扩充过的样本和 CML estimator 来估计 $β$

举个例子：
考虑一个个体被观察两次。我们首先复制 $K - 1$ 份该个体的观察值，然后每一份都在不同的截断点二分。假如一份副本 $i$ 在截断点 3 二分，那么有 $d_{i} = {1 (y_{i 1} \geq 3), 1 (y_{i 2} \geq 3)}^{'}$ ，则下一份副本 $j$ 在截断点 4 二分，有 $d_{j} = {1 (y_{i 1} \geq 4), 1 (y_{i 2} \geq 4)}^{'}$ 。

2.2.3 BUC - τ estimator

在前面的模型设定中，我们提到过每个个体的阈值 $τ_{i}$ 是不同的，且 BUC estimator 中也遵循该设定。但在学术研究中，标准的固定效应有序 Logit 模型中假设不同个体的阈值是相同的： $τ_{1 k} = τ_{2 k} = \dots = τ_{i k} = τ_{k}$ 。

因此模型设定中的 (2) 式转变为:

Pr (y_{i t} = k | x_{i t}^{'}, α_{i}) = Λ (τ_{k + 1} - x_{i t}^{'} β - α_{i}) - Λ (τ_{k} - x_{i t}^{'} β - α_{i})

如同我们在 2.1 节中所说，该公式无法单独识别出 $τ_{k}$ 和 $α_{i}$ 。因此，我们假设 $τ_{2} = 0$ 。则 $(1)$ 变为：

τ_{1} = - \infty; τ_{2} = 0; - \infty < τ_{k} < τ_{k + 1} < \infty, k = 3, \dots, K - 1; τ_{K + 1} = \infty (1)

在前面的 BUC estimator 中，每个副本内是使用同一个截断点，但在 BUC- $τ$ estimator 中，每个副本内使用随机的截断点。让 $τ_{i}^{c u t}$ 代表副本 $i$ 使用的阈值向量。则 $(3)$ 式可以重新表述为：

Pr (d_{i}^{k} | \sum_{t = 1}^{T} d_{i t} = g_{i}) = \frac{\exp {d_{i}^{^{'}} (x_{i} β - τ_{i}^{c u t})}}{\sum_{j \in B_{i}} \exp {j^{'} (x_{i} β - τ_{i}^{c u t})}} (6)

可以看到， $(6)$ 式取决于阈值。

考虑 2.2.2 中相同的例子：
假设副本 $i$ 中的第一个观察值在截断点 3 二分，第二个观察值在截断点 4 二分。那么 $d_{i} = {1 (y_{i t} \geq 3), 1 (y_{i t} \geq 4)}^{'}$ ，相应的 $τ_{i}^{c u t} = (τ_{3}, τ_{4})^{'}$ 。

于是，我们可以将上述的两个估计方法进行比较：

	BUC	BUC - τ
每个个体副本个数	$K - 1$	$(K - 1)^{T}$
样本总数	$N (K - 1)$	$N (K - 1)^{T}$

但在实际运用中， $N (K - 1)^{T}$ 超出了很多计算机的运算能力。因此，我们在程序 feologit, threshold 中，只使用一部分使用不同截断点的副本，默认为每个个体 10 个。当然，用户也可以使用 clones() 来改变数量。另外，截断点是随机挑选出来的，受随机数 seed(#) 影响。

3. 模型解释

3.1 基于方向和 compensating variation 解读 β

$β$ 的符号反映了 $x$ 的增加影响 $y$ 的累积分布变动的方向：

Pr (y_{i} t \geq k | x_{i t}, α_{i}) = Λ (x_{i t}^{'} β + α_{i} - τ_{k}) (7)

例如：

假如 $β_{l} > 0$ ，那么自变量 $x_{l}$ 的增加会减少最小截断点处的 $Pr (y_{i t} \geq 1 | x_{i t}, α_{i})$ 并增加最大截断点处的 $Pr (y_{i t} \geq K | x_{i t}, α_{i})$ 。

另外，也可以通过计算 compensating variation 来解读 $β$ ：

假如自变量 $x_{l}$ 和 $x_{r}$ 的变化使得潜在变量 $y^{*}$ 不变，进而导致因变量 $y$ 不变，那么 compensating variation 可以由相应的 $β$ 得出：增加 1 单位的 $x_{l}$ 和增加 $β_{l} / β_{r}$ 单位的 $x_{r}$ 具有相同的效果。

3.2 基于胜算比 (Odds) 解读 β

在 Logit 模型中，胜算比 (Odds) 是一个很重要的指标，它表示某一件事发生的概率与其互补事件发生的概率的比值。在本文阐述的模型中，根据 $(7)$ 式，胜算比可以表示为：

O d d s (k, x_{i t}) \equiv \frac{Pr (y_{i t} > k | x_{i t})}{Pr (y_{i t} \leq k | x_{i t})} = \exp (x_{i t}^{'} β - τ_{i k})

假设第 $l$ 个解释变量增加了 $△ x_{i t l}$ 个单位，那么胜算比的变动为：

O d d s (k, △ x_{i t l}) \equiv \frac{O d d s (k, x_{i t} + △ x_{i t l})}{O d d s (k, x_{i t})} = \exp (△ x_{i t l}^{'} β)

因此，系数 $β_{l}$ 可以解释为：

第 $x_{l}$ 每增加 1 个单位，除了最小截断点处的所有截断点处的胜算比都会变为原来的 $\exp (β_{l})$ 倍。

feologit 命令的选项 or 可以将结果呈现为 $e x p (β)$ 。例如：

. use nlswork.dta, clear

. recode hours (0/6 = 1) (7/29 = 2) (30/186 = 3), gen(hourscat)

. feologit hourscat age union msp nev_mar tenure ln_wage,or nolog group(idcode)

note: multiple positive outcomes within groups encountered.

Fixed-effects ordered logistic regression

                            N. of obs. (inc. copies) =       6874
                            N. of observations       =       6168
                            N. of panel units        =       1148
                            Wald chi2(6)             =     145.74
                            Prob > chi2              =     0.0000
Log conditional likelihood = -2398.5032  Pseudo R2   =     0.0616
                 (Std. Err. adjusted for 1,148 clusters in idcode)
------------------------------------------------------------------
         |               Robust
hourscat | Odds Ratio   Std. Err.      z    P>|z|        [95% CI]
---------+-------------------------------------------------------
     age |      0.973      0.008    -3.16   0.002   0.957   0.990
   union |      2.160      0.284     5.86   0.000   1.670   2.794
     msp |      0.472      0.069    -5.15   0.000   0.355   0.628
 nev_mar |      2.103      0.585     2.68   0.007   1.220   3.626
  tenure |      1.134      0.022     6.56   0.000   1.092   1.177
 ln_wage |      0.986      0.128    -0.11   0.912   0.764   1.271
-----------------------------------------------------------------

3.3 边际效应 (Marginal Effects) 解读

在实际应用中，我们更加关注概率值的边际效用，即一个解释变量 $x_{l}$ 的小幅度变动对 $y_{i t} = k$ 的概率值的影响，结合 $(2)$ 和 $(7)$ 式可得：

M E_{i t k l} = \frac{\partial Pr (y_{i t} = k | x_{i t}, α_{i})}{\partial x_{i t l}} = [Λ (τ_{i k + 1} - x_{i t} - α_{i}) {1 - Λ (τ_{i k + 1} - x_{i t} - α_{i})} - Λ (τ_{i k} - x_{i t} - α_{i}) {1 - Λ (τ_{i k} - x_{i t} - α_{i})}] β_{l} (8)

为了计算上式，常见方法为使用样本均值，计算出来的 Marginal Effects 称为平均边际效应 (ME at the average)：

{\bar{M E}}_{k l} = {{\bar{d}}^{k + 1} (1 - {\bar{d}}^{k + 1}) - {\bar{d}}^{k} (1 - {\bar{d}}^{k})} {\hat{β}}_{l}

其中， ${\bar{d}}^{k}$ 为 $d_{i t}^{k}$ 的样本均值。