👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：最新专题 | 计量专题 | 关于连享会
🍓 课程推荐：连享会：2025 寒假前沿班
嘉宾：杨海生，中山大学
时间：2025 年 1 月 13-24 日
咨询：王老师 18903405450（微信）
作者：李胜胜 (安徽大学)
邮箱：lishengsheng2016@126.com
编者按：本文摘译自下文，特此致谢！
Source：Correia S, Guimarães P, Zylkin T. Fast Poisson estimation with high-dimensional fixed effects[J]. The Stata Journal, 2020, 20(1): 95-115. -PDF-
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
目录
1. 引言
2. 估计方法
2.1 IRLS 算法
2.2 泊松回归
3. ppmlhdfe 命令
3.1 命令安装：
3.2 语法结构：
4. Stata Stata 范例
4.1 范例 1：模拟数据
4.2 范例 2：Stata 手册例子
4.3 范例 3：贸易引力模型
4.4 范例 4：结果输出及展示
5. 参考文献
6. 相关推文
1. 引言
泊松回归现已成为「计数数据模型」的标准方法。与普通最小二乘法一样，泊松回归的一致性估计也需要正确说明因变量的条件均值假设 (Gourieroux 等，1984)。在该条件下，泊松回归就能转化为泊松伪最大似然 (PPML) 回归。Gourieroux 等 (1984) 通过放松因变量分布的假设，使泊松回归不再局限于计数数据。也因此，泊松伪最大似然 (PPML) 回归可以应用于任何非负的因变量。
实际上，当非负数据中存在较多 0 的时候，PPML 似乎是最好的选择。并且，这种情况常发生于以下几个领域，例如公司研发支出、专利引用计数、日常产品销售、医生就诊次数、公司信贷量、拍卖竞拍者数量和跨地区通勤者数量。
但是，学者更多使用线性对数回归，即使 PPML 更合理。一种可能的解释是，学者可以轻松估计控制多个固定效应的线性回归。尤其是伴随着大型面板数据集的可用性不断提高，以及具有高维固定效应 (HDFE) 的线性回归模型估计技术的进步，使得学者能够控制更多异质性来源。
需要指出的是，PPML 也可以轻松的实现 HDFE。为此，本文将介绍新命令 ppmlhdfe。相对于用于现有 HDFE 非线性估计算法，ppmlhdfe 消除了一些不必要的步骤，提高了参数的计算速度。
2. 估计方法
2.1 IRLS 算法
GLM 是基于 Nelder 和 Wedderburn (1972) 引入的指数分布族的一类回归模型，主要包括流行的非线性回归模型，例如 Logit、Probit、Cloglog 和 Poisson。Hardin 和 Hilbe (2018) 将指数族定义如下：
其中，、 和 是特定函数， 和 是参数。对于这些模型而言：
给定一组具有 个独立的观察值，每个观察值均由 索引，可以通过链接函数 将期望值与一组协变量 () 关联。更具体地说，假设下式：
则 GLM 的似然函数可以写为：
通过求解 (伪) 似然最大化的一阶条件，得到 的估计。应用高斯-牛顿算法和期望 Hessian 函数导出修正方程为：
其中 是解释变量矩阵， 是权重矩阵， 是因变量的变换， 是迭代索引。式 (1) 表明 的估计值是通过加权最小二乘法递归而获得，这种方法称之为 IRLS。
2.2 泊松回归
泊松回归定义如下：
实现 IRLS 的回归权重可以简化为：
而中间回归的因变量则为:
因此，在 (1) 中实现 IRLS 更新回归仅需要计算在先前迭代中获得拟合值的向量即可。
处理 HDFE
在 HDFE 下，IRLS 的困难是 可能包含很多固定效应，导致无法直接计算 。解决方案是使用一个替代的公式，即只估计非固定效应协变量的系数 (例如，)，从而降低维数。由于式 (1) 是加权线性回归，因此可以依靠 FWL 定理来探讨固定效应。这意味着代替式 (1)，可以使用下式：
其中， 和 分别是主协变量矩阵 和因变量 的变换后的加权。此外，FWL 定理还暗示从式 (4)计算的残差与从式 (1) 计算的残差相同。这意味着可以使用以下方法对 和 执行所需的更新：
其中， 是一个向量，该向量收集使用式 (4) 计算的残差。然后，与原始 IRLS 循环一样， 和 的新值直接从式 (2) 和 (3) 开始。其次，这也意味着，一旦 收敛到正确的估计值 ，式 (4) 中加权最小二乘回归的估计方差-协方差矩阵将成为 正确的方差-协方差矩阵。
加速 HDFE–GIRLS
命令 poi2hdfe 可以实现式 (4)，同时使用 reghdfe 作为运行 HDFE 加权最小二乘回归。这是一个密集型计算过程，需要在每次 IRLS 迭代中估算 HDFE 回归模型。但是，ppmlhdfe 中有多种变通办法，可以使其效率更高。例如，ppmlhdfe 直接嵌入 reghdfe 的 Mata 中。因此利用了这样一个事实，它们在每个 IRLS 迭代中保持不变，某些计算只需要执行一次。但最显著的速度改进来自于对标准 HDFE–IRLS 算法的改进，该算法旨在减少对 reghdfe 的调用次数。
最大似然估计的存在性
Santos Silva 和 Tenreyro (2010，2011) 指出，对于某些数据配置，可能不存在泊松回归的最大似然估计 (MLE)，进而导致估计可能无法收敛或收敛到错误的值。在泊松回归的情况下，如果对数似然随着一个或多个系数趋于无穷大而单调增加，则会发生这种情况。Santos Silva 和 Tenreyro (2010) 认为发生这种情况的主要原因是变量间的多重共线性。为此，他们建议排除有问题的变量。但是，排除哪个回归变量是一个模棱两可的决定，可能会影响其余参数的识别。此外，在具有多个 HDFE 的泊松模型中，该策略甚至不可行。
Correia 等 2019 讨论了各种 GLM 模型估计中的必要条件和充分条件，并表明在泊松回归情况下，如果从样本中删除一些观察值，可能得到总 MLE 估计。这些单独的观测值不传递估计过程的相关信息，可以安全地丢弃。同时，删除这些观察值后，某些回归变量将产生共线性，因此也必须删除。此外，Correia 等 2019 提出了一种识别分离观察结果的方法，并且即使在 HDFE 存在情况下，也可以成功运行。
3. ppmlhdfe 命令
3.1 命令安装：
ssc install ppmlhdfe, replace
3.2 语法结构：
ppmlhdfe depvar [indepvars] [if] [in] [weight] , ///
[absorb(absvars)] [options]
depvar：被解释变量；indepvars：解释变量；absorb(absvars)：要吸收的分类变量 (固定效应)，也允许单独的斜率；absorb(..., savefe)：使用 hdfe # 保存所有固定效应估计值。[options]选项：exposure(varname)：在系数约束为 1 的模型中包含 ln(varname)；offset(varname)：在系数约束为 1 的模型中包含 varname；d(newvar)：将固定效应之和另存为newvar；d：如上，但变量将另存为 _ppmlhdfe_d；vce(vcetype)：vcetype 可以是 robust 或聚类；verbose(#)：显示调试信息量；nolog：隐藏迭代日志；tolerance(#)：收敛标准，默认为tolerance(1e-8)；guess(string)：设置用于设置初始值的规则；eform：报告指数系数；irr：eform 的同义词；separation(string)：用于删除分离的观测值及其相关回归变量的算法；maxiteration(#)：指定最大迭代次数；keepsingletons：不要删除单例组；version：报告ppmlhdfe的版本号和日期；display_options：控制回归表的选项，如置信水平、数字格式等。
4. Stata Stata 范例
4.1 范例 1：模拟数据
从一个简单的例子开始，我们来说明 ppmlhdfe 方法处理 MLE 估计不存在时的优势。如前所述，ppmlhdfe 注意识别分离的观测值，然后以保证存在有意义的 MLE 方式限制样本。
以下数据包括六个观察值和三个解释变量：
input y x1 x2 x3
0 1 2 1
0 0 0 2
0 2 3 3
1 1 2 4
2 2 4 5
3 1 2 6
end
如果试图用 glm 命令估计泊松回归，Stata 就无法收敛。poisson 命令可以得到三个解释变量系数的估计值，但是会发现估计结果并不可靠。
. poisson y x1 x2 x3, nolog
Poisson regression Number of obs = 6
LR chi2(3) = 8.89
Prob > chi2 = 0.0308
Log likelihood = -4.0415302 Pseudo R2 = 0.5237
-----------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------+----------------------------------------------------------------
x1 | -31.29521 8467.059 -0.00 0.997 -16626.43 16563.83
x2 | 15.84339 4233.529 0.00 0.997 -8281.722 8313.408
x3 | .7970409 .4608654 1.73 0.084 -.1062388 1.70032
_cons | -4.032453 2.868066 -1.41 0.160 -9.653759 1.588853
-----------------------------------------------------------------------
使用命令 ppml 识别存在问题的数据，并删除变量 。但是，拟合回归仍然有问题。
. ppml y x1 x2 x3
note: checking the existence of the estimates
note: starting ppml estimation
Iteration 1: deviance = 5.984675
Iteration 2: deviance = 5.625005
Iteration 3: deviance = 5.538839
Iteration 4: deviance = 5.507367
Iteration 5: deviance = 5.49579
Iteration 6: deviance = 5.49153
Iteration 7: deviance = 5.489963
Iteration 8: deviance = 5.489387
Iteration 9: deviance = 5.489175
Iteration 10: deviance = 5.489097
Iteration 11: deviance = 5.489068
Iteration 12: deviance = 5.489058
Iteration 13: deviance = 5.489054
Iteration 14: deviance = 5.489052
Iteration 15: deviance = 5.489052
Warning: variance matrix is nonsymmetric or highly singular
Number of parameters: 3
Number of observations: 6
Number of observations dropped: 0
Pseudo log-likelihood: -6.5473014
R-squared: .3399843
WARNING: The model appears to overfit some observations with y=0
-----------------------------------------------------------------------
| Robust
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------+----------------------------------------------------------------
x1 | -32.31708 . . . . .
x2 | 16.58591 . . . . .
_cons | -.8167293 . . . . .
-----------------------------------------------------------------------
Number of regressors dropped to ensure that the estimates exist: 1
Dropped variables: x3
Option strict is off
如果使用 strict 选项运行，则 ppml 会删除所有回归变量。而 ppmlhdfe 方法不同，与其寻找有问题的回归变量，不如寻找有问题的观察值。在这种情况下，ppmlhdfe 删除了第三个观察值，然后丢弃 以避免完全多重共线，此时估计结果会更合理。
. ppmlhdfe y x1 x2 x3, nolog
(simplex method dropped 1 separated observation)
note: 1 variable omitted because of collinearity: x2
Converged in 6 iterations and 6 HDFE sub-iterations (tol = 1.0e-08)
PPML regression No. of obs = 5
Residual df = 2
Wald chi2(2) = 50.78
Deviance = .4775093816 Prob > chi2 = 0.0000
Log pseudolikelihood = -4.041530113 Pseudo R2 = 0.4532
-----------------------------------------------------------------------
| Robust
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------+----------------------------------------------------------------
x1 | .3914642 .1733025 2.26 0.024 .0517976 .7311309
x2 | 0 (omitted)
x3 | .7969293 .1582404 5.04 0.000 .4867839 1.107075
_cons | -4.031679 1.119577 -3.60 0.000 -6.226011 -1.837348
-----------------------------------------------------------------------
4.2 范例 2：Stata 手册例子
接下来，我们使用固定效应选项的命令 xtpoisson 来复制 Stata 16 手册中的示例 1。本示例使用的数据为 ships.dta，并以此来估计多个回归变量上船舶事故数量的泊松回归。它将可变船舶作为固定效应，即控制五种类型的船舶。通过控制 exposure(service) 估计回归，并将系数报告为发生率。
. webuse ships, clear
. xtpoisson acc op_75_79 co_65_69 co_70_74 co_75_79, ///
exposure(service) irr fe nolog
Conditional fixed-effects Poisson regression Number of obs = 34
Group variable: ship Number of groups = 5
Obs per group:
min = 6
avg = 6.8
max = 7
Wald chi2(4) = 48.44
Log likelihood = -54.641859 Prob > chi2 = 0.0000
--------------------------------------------------------------------------
accident | IRR Std. Err. z P>|z| [95% Conf. Interval]
-------------+------------------------------------------------------------
op_75_79 | 1.468831 .1737218 3.25 0.001 1.164926 1.852019
co_65_69 | 2.008002 .3004803 4.66 0.000 1.497577 2.692398
co_70_74 | 2.26693 .384865 4.82 0.000 1.625274 3.161912
co_75_79 | 1.573695 .3669393 1.94 0.052 .9964273 2.485397
ln(service) | 1 (exposure)
--------------------------------------------------------------------------
用 ppmlhdfe 可以获得等效的结果，语法为：
. ppmlhdfe acc op_75_79 co_65_69 co_70_74 co_75_79, ///
absorb(ship) exposure(service) irr nolog
Converged in 6 iterations and 6 HDFE sub-iterations (tol = 1.0e-08)
HDFE PPML regression No. of obs = 34
Absorbing 1 HDFE group Residual df = 25
Wald chi2(4) = 111.06
Deviance = 38.69505154 Prob > chi2 = 0.0000
Log pseudolikelihood = -68.28077143 Pseudo R2 = 0.8083
-----------------------------------------------------------------------
| Robust
accident | exp(b) Std. Err. z P>|z| [95% CI]
-------------+---------------------------------------------------------
op_75_79 | 1.468831 .1484359 3.80 0.000 1.204902 1.790572
co_65_69 | 2.008002 .2202475 6.36 0.000 1.619572 2.489592
co_70_74 | 2.26693 .3256501 5.70 0.000 1.710649 3.004107
co_75_79 | 1.573695 .3117262 2.29 0.022 1.067358 2.320232
_cons | .0011254 .0001061 -72.03 0.000 .0009356 .0013538
ln(service) | 1 (exposure)
-----------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
Absorbed FE | Categories - Redundant = Num. Coefs |
-------------+---------------------------------------|
ship | 5 0 5 |
-----------------------------------------------------+
在这里吸收船舶作为固定效应，有几点值得一提，如预期的那样，变量的估计系数与使用 xtpoisson 命令获得的系数相同。但是，标准误的估计结果有所不同，因为默认情况下，ppmlhdfe 报告的是可靠的标准误。这两个命令显示的对数似然值也不同，命令 xtpoisson 报告条件对数似然值，而 ppmlhdfe 报告实际的泊松对数似然值。
考虑到我们使用的是一个小数据集，可以使用 poisson 命令复制ppmlhdfe 获得的结果，如下所示：
. ppmlhdfe acc op_75_79 co_65_69, ///
absorb(ship co_70_74 co_75_79) ///
exposure(service) irr nolog
值得注意的是，ppmlhdfe 可以吸收任何类别变量作为固定效应。例如，如果仅对 op_75_79 和 co_65_69 的系数感兴趣，则可以将 ship、co_70_74 和 co_ 75_79 吸收。
. ppmlhdfe acc op_75_79 co_65_69, ///
absorb(ship co_70_74 co_75_79) ///
exposure(service) irr nolog
Converged in 7 iterations and 23 HDFE sub-iterations (tol = 1.0e-08)
HDFE PPML regression No. of obs = 34
Absorbing 3 HDFE groups Residual df = 25
Wald chi2(2) = 71.60
Deviance = 38.69505154 Prob > chi2 = 0.0000
Log pseudolikelihood = -68.28077143 Pseudo R2 = 0.8083
-------------------------------------------------------------------------
| Robust
accident | exp(b) Std. Err. z P>|z| [95% CI]
-------------+-----------------------------------------------------------
op_75_79 | 1.468831 .1484359 3.80 0.000 1.204902 1.790572
co_65_69 | 2.008002 .2202475 6.36 0.000 1.619572 2.489592
_cons | .0015435 .0001204 -83.00 0.000 .0013247 .0017984
ln(service) | 1 (exposure)
-------------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
Absorbed FE | Categories - Redundant = Num. Coefs |
-------------+---------------------------------------|
ship | 5 0 5 |
co_70_74 | 2 1 1 |
co_75_79 | 2 1 1 ?|
-----------------------------------------------------+
? = number of redundant parameters may be higher
4.3 范例 3：贸易引力模型
使用 ppml_ panel_sg 命令提供的辅助数据和示例来拟合引力模型。该数据集包含 35 个国家 1986 年至 2004 年的年度双边贸易数据。目的是估计自由贸易协定变量 fta 对贸易的影响。在本例中，我们希望控制国家对固定效应 (country-pair fixed effects) 和国家时间固定效应 (对于进口国和出口国)。此外，希望标准误聚类在国家对的级别上。
* 数据下载地址：
* https://gitee.com/arlionn/data/blob/master/data01/EXAMPLE_TRADE_FTA.dta
use http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA if category=="TOTAL", clear
egen imp=group(isoimp)
egen exp=group(isoexp)
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) cluster(imp#exp) nolog
Converged in 11 iterations and 35 HDFE sub-iterations (tol = 1.0e-08)
HDFE PPML regression No. of obs = 5,950
Absorbing 3 HDFE groups Residual df = 1,189
Statistics robust to heteroskedasticity Wald chi2(1) = 21.04
Deviance = 377332502.3 Prob > chi2 = 0.0000
Log pseudolikelihood = -188710931.7 Pseudo R2 = 0.9938
Number of clusters (imp#exp)= 1,190
(Std. Err. adjusted for 1,190 clusters in imp#exp)
-------------------------------------------------------------------
| Robust
trade | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------+-----------------------------------------------------------
fta | .1924455 .0419527 4.59 0.000 .1102197 .2746713
_cons | 16.45706 .0217308 757.32 0.000 16.41447 16.49965
-------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
Absorbed FE | Categories - Redundant = Num. Coefs |
-------------+---------------------------------------|
imp#year | 175 0 175 |
exp#year | 175 5 170 |
imp#exp | 1190 1190 0 *|
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation
因为交易示例包括 HDFE，所以它也为我们提供了一个机会，可以从编程的角度解开 ppmlhdfe 的工作原理和背后运行机制。ppmlhdfe 本身是一个复杂的命令，主要在 Mata 中进行编码，并尽可能利用 reghdfe 现有内部工作原理。但是，用于实现 HDFE–IRLS 基本算法具有很高的可移植性，并且可以使用任何已经为内部转换和标准加权回归提供可靠算法的语言进行编程。
4.4 范例 4：结果输出及展示
本示例将展示如何将 ppmlhdfe 与 esttab 命令结合来构建发布高质量的回归表格。为了构造固定效应的标签，使用了 reghdfe 附带的 estfe 命令。该示例由一个包含 1,000,000 多个观察值的数据组成，其中包含 1991 年至 2009 年间收集的所有 SSC 引文。该数据具有期刊 ID、文章标题、作者人数、两个主要的 JEL 分类代码、出版年份、文章类型、以及每年收集的引用次数。但是，所有观测值中有将近 58％ 为零。因变量是引文的数量，我们主要的关注点在于一旦控制了各种可以控制的变量，一篇文章中作者数量是否真的提升了引文数量。
* 数据下载地址：
* https://gitee.com/arlionn/data/blob/master/data01/citations_example.rar
use citations_example, clear
estimates clear
ppmlhdfe cit nbaut, absorb(issn type jel2 pubyear)
eststo
ppmlhdfe cit nbaut, absorb(issn#c.year type jel2 pubyear)
eststo
ppmlhdfe cit nbaut, absorb(issn#year type jel2 pubyear)
eststo
estfe *, labels(type "Article type FE" jel2 "JEL code FE" pubyear "Publication year FE" ///
issn "ISSN FE" issn#c.year "Year trend by ISSN" issn#year "ISSN-Year FE")
esttab, indicate(`r(indicate_fe)', ///
labels("Yes" "")) b(3) se(3) ///
varwidth(25) label compress ///
stat(N ll, fmt(%9.0fc %9.2fc)) ///
se starlevels(* 0.1 ** 0.05 *** 0.01)
输出结果如下：
----------------------------------------------------------------
(1) (2) (3)
cit cit cit
----------------------------------------------------------------
nbaut 0.190*** 0.190*** 0.189***
(0.003) (0.003) (0.003)
Constant 0.054*** 0.081*** 0.110***
(0.006) (0.006) (0.006)
Article type FE Yes Yes Yes
JEL code FE Yes Yes Yes
Publication year FE Yes Yes Yes
ISSN FE Yes
Year trend by ISSN Yes
ISSN-Year FE Yes
----------------------------------------------------------------
N 1,083,701 1,083,701 1,080,051
ll -1.71e+06 -1.69e+06 -1.66e+06
----------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01
结果表明作者数量对被引次数具有显著正影响。
5. 参考文献
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
Gourieroux C, Monfort A, Trognon A. Pseudo maximum likelihood methods: Theory[J]. Econometrica: journal of the Econometric Society, 1984: 681-700. -PDF- Correia S, Guimarães P, Zylkin T. Fast Poisson estimation with high-dimensional fixed effects[J]. The Stata Journal, 2020, 20(1): 95-115. -PDF-
6. 相关推文
Note：产生如下推文列表的命令为：
lianxh 固定效应 泊松
安装最新版lianxh命令：
ssc install lianxh, replace
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
专题：面板数据 ocmt：高维固定效应模型的变量筛选问题 Stata新命令：ppmlhdfe-面板计数模型-多维固定效应泊松估计 Stata：非对称固定效应模型 reghdfe：多维面板固定效应估计 专题：倍分法DID Stata：双重差分的固定效应模型-(DID) 专题：内生性-因果推断 用FE-固定效应模型能做因果推断吗？ 专题：Probit-Logit feologit：固定效应有序Logit模型 ZIP-too many Zero：零膨胀泊松回归模型
尊敬的老师 / 亲爱的同学们：
连享会致力于不断优化和丰富课程内容，以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求，我们诚挚地邀请您参与到我们的课程规划中来。 请您在下面的问卷中，分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源，将直接影响到我们课程的改进和创新。 我们期待您的反馈，因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间，与我们共同塑造更加精彩的学习旅程！https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见！
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。
