温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装命令如下:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
⛳ Stata 系列推文:
作者:李胜胜 (安徽大学)
邮箱: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-
目录
泊松回归现已成为「计数数据模型」的标准方法。与普通最小二乘法一样,泊松回归的一致性估计也需要正确说明因变量的条件均值假设 (Gourieroux 等,1984)。在该条件下,泊松回归就能转化为泊松伪最大似然 (PPML) 回归。Gourieroux 等 (1984) 通过放松因变量分布的假设,使泊松回归不再局限于计数数据。也因此,泊松伪最大似然 (PPML) 回归可以应用于任何非负的因变量。
实际上,当非负数据中存在较多 0 的时候,PPML 似乎是最好的选择。并且,这种情况常发生于以下几个领域,例如公司研发支出、专利引用计数、日常产品销售、医生就诊次数、公司信贷量、拍卖竞拍者数量和跨地区通勤者数量。
但是,学者更多使用线性对数回归,即使 PPML 更合理。一种可能的解释是,学者可以轻松估计控制多个固定效应的线性回归。尤其是伴随着大型面板数据集的可用性不断提高,以及具有高维固定效应 (HDFE) 的线性回归模型估计技术的进步,使得学者能够控制更多异质性来源。
需要指出的是,PPML 也可以轻松的实现 HDFE。为此,本文将介绍新命令 ppmlhdfe
。相对于用于现有 HDFE 非线性估计算法,ppmlhdfe
消除了一些不必要的步骤,提高了参数的计算速度。
GLM 是基于 Nelder 和 Wedderburn (1972) 引入的指数分布族的一类回归模型,主要包括流行的非线性回归模型,例如 Logit、Probit、Cloglog 和 Poisson。Hardin 和 Hilbe (2018) 将指数族定义如下:
其中,
给定一组具有
则 GLM 的似然函数可以写为:
通过求解 (伪) 似然最大化的一阶条件,得到
其中
泊松回归定义如下:
实现 IRLS 的回归权重可以简化为:
而中间回归的因变量则为:
因此,在 (1) 中实现 IRLS 更新回归仅需要计算在先前迭代中获得拟合值
处理 HDFE
在 HDFE 下,IRLS 的困难是
其中,
其中,
加速 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 存在情况下,也可以成功运行。
ssc install ppmlhdfe, replace
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
:控制回归表的选项,如置信水平、数字格式等。
从一个简单的例子开始,我们来说明 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
-----------------------------------------------------------------------
接下来,我们使用固定效应选项的命令 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
使用 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 基本算法具有很高的可移植性,并且可以使用任何已经为内部转换和标准加权回归提供可靠算法的语言进行编程。
本示例将展示如何将 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
结果表明作者数量对被引次数具有显著正影响。
Note:产生如下推文列表的命令为:
lianxh 固定效应 泊松
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh