Stata：一文读懂两部模型-twopm
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会课程 · 基于机器学习的因果推断方法
作者：王恩泽 (武汉大学经济与管理学院)
邮箱：enzewang2008@163.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
编者按：本文主要摘译自下文，特此致谢！
Source：Belotti F, Deb P, Manning W G, et al. twopm: Two-part models[J]. The Stata Journal, 2015, 15(1): 3-20. -PDF-
目录
1. 简介
1.1 什么是两部分模型？
1.2 twopm 命令简介
2. twopm 命令的语法格式
2.1 语法格式 1 (两部分为相同的回归模型)
2.2 语法格式 2 (两部分为不同的回归模型)
2.3 语法选项详解
2.4 使用 twopm 语法后的估计与预测
3. Stata 范例
3.1 范例 1
3.2 范例 2
3.3 总结
4. 结语
5. 相关推文
1. 简介
1.1 什么是两部分模型？
在实证分析的过程中，我们经常会遇到混合型 (离散-连续) 随机变量，这类变量往往有如下基本特征：(1) ≥ 0，(2) = 0。换句话说，该类变量往往由大量的 值与大于 的值组成。由于大量 值的存在，单一的指数模型往往无法对这类数据进行较好地拟合。为此，我们需要借助于两部分模型 (Two-part models) 对该类数据进行处理。
具体而言，两部分模型借助一个二元选择模型获得观测值为零值的概率。随后，以观测到正值为条件。一个合适的回归模型被应用于拟合正值。在本篇推文中，我们介绍了 twopm 命令，该命令可被应用于估计两部分模型。同时，借助该命令，我们也能够计算使用两部分模型得出的预测值与边际效应。
以下是 ChatGPT 关于两部分模型 (Two-part models) 的介绍：
Two-part models（两部分模型）是一种统计模型，由两个部分组成，用于描述具有大量零值或较小计数的数据。这种模型的两个部分通常是一个二项分布或泊松分布的概率质量函数和一个正态分布的概率密度函数。
Two-part models通常用于医疗保健和社会科学领域的数据分析，例如用于描述医疗服务利用、药物使用等问题。这种模型最常见的用途是对一些稀有事件（如住院、手术等）进行建模。
Two-part models中的第一部分通常描述了事件的发生概率，因为大部分数据都是零。第二部分描述了对于发生了事件的个体，其数值数据如何分布。这种模型通常使用最大似然估计进行参数估计。
Two-part models具有许多优点，例如可以对具有复杂数据结构的数据进行建模，可以更好地描述实际数据的特征。但是，它们也有一些缺点，例如可能需要更多的计算资源和更复杂的模型选择过程。
以下是一些关于Two-part models的例子：
医疗保健利用模型：Two-part models常用于医疗保健领域的数据分析，例如描述患者住院率、手术率、就诊次数等问题。该模型的第一部分可以描述患者发生住院或手术的概率，第二部分可以描述住院或手术的费用分布。
药物使用模型：Two-part models也可用于描述药物使用情况，例如描述患者是否使用某种药物和他们使用药物的频率和剂量。该模型的第一部分可以描述患者是否使用药物的概率，第二部分可以描述药物剂量的分布。
家庭收入模型：Two-part models还可用于社会科学领域的数据分析，例如描述家庭收入情况。该模型的第一部分可以描述家庭收入为零的概率，第二部分可以描述家庭收入非零的分布。
保险索赔模型：Two-part models也可用于描述保险索赔的情况，例如描述保险索赔是否发生以及发生时索赔的金额。该模型的第一部分可以描述索赔是否发生的概率，第二部分可以描述索赔金额的分布。
这些例子仅是Two-part models应用的一小部分，该模型可用于任何具有大量零值或较小计数的数据集。
1.2 twopm 命令简介
twopm 命令允许我们在第一步中使用 logit 或 probit 模型对数据进行拟合，在第二步中使用 OLS 或者 GLM 回归对数据进行拟合。相比于分别对两部分模型进行估计，twopm 命令具有如下优势：
第一，它可以适应复杂的调查设计，进而对参数估计和这些估计的标准误差进行调整。复杂的调查设计在大型调查中很常见，然而，忽略调查结构可能导致总体参数的估计偏差。 第二，该命令可以对两部分模型中每一步都存在变量的联合显著性进行检验。 第三，使用该命令后，借助 predict和margins可以方便的计算被解释变量的预测值与解释变量的边际效应。值得注意的是，此处的预测值是针对全样本而言的，而非仅仅针对第二步预测过程中包含的正数观测值。第四，该命令可以对估计尺度进行再转换。具体而言，在第二步估计时若使用 OLS 对 进行回归，该程序能够提供以 尺度 ( 原始尺度 ) 为基础预测的估计值。 第五，该命令会在考虑模型的两个部分、复杂的调查设计、及基于增量方法 (delta method) 的稳健标准误的基础上自动计算预测值和边际效应的标准误差。
2. twopm 命令的语法格式
该命令具有两种语法格式，第一种语法格式自动在第一和第二步中定义相同的回归模型和函数形式；第二种语法格式允许我们自主设定第一和第二步中所使用的的回归模型。
ssc install twopm , replace
2.1 语法格式 1 (两部分为相同的回归模型)
twopm depvar [indepvars] [if] [in] [weight],
firstpart(f_options) secondpart(s_options)
[vce(vcetype) robust cluster(cluster)
suest level(#) nocnsreport display_options]
2.2 语法格式 2 (两部分为不同的回归模型)
twopm equation1 equation2 [indepvars] [if] [in] [weight],
firstpart(f_options) secondpart(s_options)
[vce(vcetype) robust cluster(cluster)
suest level(#) nocnsreport display_options]
其中，equation1 和 equation2 被定义为如下形式：
(depvar [=] [indepvars])
2.3 语法选项详解
firstpart(f_options)：设定模型第一步估计采用的形式。(f_options)处填写logit或probit。secondpart(s_options)：设定模型第二步估计采用的形式。(s_options)处填写regress或glm。vce(vcetype)：设定模型所报告标准误的类型。五种类型可选：包括从渐近理论推导出的类型、对某些类型的错误识别具有稳健性的类型、允许组内相关的类型、以及使用 bootstrap 或 jackknife 方法的类型。默认类型为vce(conventional)，即在模型的第一步和第二步中使用常规的方差估计。robust：相当于使用vce(robust)。cluster (clustvar)：相当于使用vce(cluster clustvar)。suest：结合了模型第一部分和第二部分的估计结果，该选项往往填写test或testnl进而用来检验某一变量的联合显著性。
2.4 使用 twopm 语法后的估计与预测
predict [type] newvar [if] [in], [{normal|duan} scores nooffset]
该语法被应用于计算预测值或者 E(y|x) 的估计值。
predict [type] {stub*|newvar1 ... newvarq} [if] [in], scores
该语法被应用于计算等式得分 (equation-level scores)。
在上述两种语法格式中，若在 [if] 处填写 if e(sample)，则意味着预测值与计算所得的等式得分都被限制于估计的样本中。在使用对数形式的被解释变量进行两部分回归后，如下选项将可以使用：
normal：使用正态理论重新转换进而获得拟合值。duan：使用 Duan (1983) 的 smearing 再转换进而获取拟合值。scores：为模型的每个部分创建得分变量。由于模型第二部分的得分只有在估计 分样本时有意义，故计算第二部分的得分变量时样本被自动限于估计的分样本。nooffset：设定计算预测值或得分时忽略设定模型时定义的偏移变量或暴露变量。
3. Stata 范例
在本部分，我们使用 2004 年医疗支出调查数据，展示了两部分模型的两个例子。具体而言，我们使用两部分模型的两种常见设定形式来估计医疗花费总支出的预测值，并计算年龄和性别对其的边际影响。twopm 命令与复杂的调查设计命令兼容，因此在读入数据后，我们使用 svyset 命令将数据设定为问卷调查的形式。
首先，导入数据并分析数据的基本统计特征，命令如下：
. lxhuse meps_ashe_subset5, clear
. svyset [pweight=wtdper], strata(varstr) psu(varpsu)
. svy: mean exp_tot age female
Survey: Mean estimation
Number of strata = 203 Number of obs = 19,386
Number of PSUs = 448 Population size = 187,973,715
Design df = 245
--------------------------------------------------------------
| Linearized
| Mean std. err. [95% conf. interval]
-------------+------------------------------------------------
exp_tot | 3838.939 99.945 3642.078 4035.801
age | 45.791 0.229 45.339 46.243
female | 0.520 0.003 0.514 0.526
--------------------------------------------------------------
3.1 范例 1
在本例中，我们在两部分模型的第一步中使用了 probit 模型，在第二步中使用了对数链接 (log link) 和伽马分布 (gamma distribution) 的 GLM 估计。
. *-Two-part model, with probit first part and GLM second part
. svy: twopm exp_tot c.age i.female, ///
> firstpart(probit) ///
> secondpart(glm, family(gamma) link(log))
Survey data analysis
Number of strata = 203 Number of obs = 19,386
Number of PSUs = 448 Population size = 187,973,715
Design df = 245
F(2, 244) = 671.26
Prob > F = 0.0000
------------------------------------------------------------------------------
| Linearized
exp_tot | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
probit |
age | 0.025 0.001 31.65 0.000 0.024 0.027
1.female | 0.564 0.027 20.76 0.000 0.511 0.618
_cons | -0.239 0.039 -6.12 0.000 -0.315 -0.162
-------------+----------------------------------------------------------------
glm |
age | 0.029 0.001 22.19 0.000 0.026 0.031
1.female | 0.200 0.054 3.70 0.000 0.093 0.306
_cons | 6.804 0.087 78.65 0.000 6.633 6.974
------------------------------------------------------------------------------
在上述命令中，exp_tot 为被解释变量医疗支出，c.age 与 i.female 分别表示关注变量年龄与性别。firstpart(probit) 表示在第一步中选用 probit 模型进行数据估计，secondpart(glm, family(gamma) link(log)) 表示在第二步中选用对数链接 (log link) 和伽马分布 (gamma distribution) 的 GLM 模型对数据进行拟合。
由结果可知，age 和 female 变量的系数均在 1% 的水平上显著为正。一方面，观察 probit 模型所示结果可知，医疗支出的概率随着年龄的增加而增长；同时，与男性相比，女性有更大的可能性花费至少 1 美元医疗支出。另一方面，观察 glm 模型所示结果可知，在任何医疗支出水平上，医疗支出的数量均随着年龄的增加而增长；同时，与男性相比，在任何支出水平上，女性更可能比男性支出的更多。
随后，我们可以使用 margins 命令进而对总花费支出进行预测。
. margins
Predictive margins
Number of strata = 203 Number of obs = 19,386
Number of PSUs = 448 Population size = 187,973,715
Model VCE: Linearized Design df = 245
Expression: twopm combined expected values, predict()
------------------------------------------------------------------------------
| Delta-method
| Margin std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | 3870.714 94.987 40.75 0.000 3683.619 4057.809
------------------------------------------------------------------------------
由结果可知，预测的总医疗支出约为 3871 美元每年，该数值与实际的年度医疗支出均值 3839 美元较为接近。接下来，借助 margins 命令，我们可以在两部分模型总体的基础上对不同变量的总体边际效应进行计算。
. margins, dydx(*)
Average marginal effects
Number of strata = 203 Number of obs = 19,386
Number of PSUs = 448 Population size = 187,973,715
Model VCE: Linearized Design df = 245
Expression: twopm combined expected values, predict()
dy/dx wrt: age 1.female
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | 127.832 6.373 20.06 0.000 115.280 140.385
1.female | 1139.541 186.779 6.10 0.000 771.642 1507.439
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
上表表明，年龄的边际效应约为每年 128 美元，即年龄增加一岁，每年的医疗支出增加约 128 美元；与此同时，女性医疗支出超出男性 1140 美元。
鉴于边际效应可能会随着年龄的变化产生变化，我们接下来分别计算四个年龄 (20，40，60，80) 处的边际效应。例如，尽管女性在所有的年龄段都比男性医疗支出更多，但年长女性的这种差异比年轻女性大得多。
. margins, dydx(*) at(age=(20(20)80))
Conditional marginal effects
Number of strata = 203 Number of obs = 19,386
Number of PSUs = 448 Population size = 187,973,715
Model VCE: Linearized Design df = 245
Expression: twopm combined expected values, predict()
dy/dx wrt: age 1.female
1._at: age = 20
2._at: age = 40
3._at: age = 60
4._at: age = 80
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age |
_at |
1 | 51.359 1.358 37.83 0.000 48.685 54.032
2 | 95.643 3.141 30.45 0.000 89.457 101.829
3 | 169.311 9.603 17.63 0.000 150.396 188.226
4 | 295.702 24.667 11.99 0.000 247.115 344.288
-------------+----------------------------------------------------------------
0.female | (base outcome)
-------------+----------------------------------------------------------------
1.female |
_at |
1 | 589.644 60.456 9.75 0.000 470.564 708.723
2 | 942.570 127.344 7.40 0.000 691.741 1193.398
3 | 1431.437 260.978 5.48 0.000 917.389 1945.484
4 | 2228.771 505.608 4.41 0.000 1232.878 3224.664
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
有时，我们也想了解是否某个解释变量在两部分模型的两部分是联合显著的。在本例中，年龄和性别在模型的每一个部分都是显著的，所以他们分别在两部分中联合显著并不令人吃惊。
. test age
Adjusted Wald test
( 1) [probit]age = 0
( 2) [glm]age = 0
F( 2, 244) = 803.99
Prob > F = 0.0000
. test 1.female
Adjusted Wald test
( 1) [probit]1.female = 0
( 2) [glm]1.female = 0
F( 2, 244) = 226.39
Prob > F = 0.0000
3.2 范例 2
在本例中，我们在两部分模型的第一步中使用了 logit 模型，在第二步中对对数形式的 使用了 OLS 估计。对于将数据向原始尺度的再转换而言，我们不要求对数尺度误差项具有正态分布的限制性假设。因为这一假设通常会导致估计的条件均值与边际效应有更大的误差。因此，我们使用了 Duan (1983) 的 smearing 估计，twopm 命令将在使用后延命令时自动使用 smearing 估计。
. twopm exp_tot c.age i.female, firstpart(logit) secondpart(regress, log)
Fitting OLS regression for second part:
Two-part model
------------------------------------------------------------------------------
Log pseudolikelihood = -37216.38 Number of obs = 19386
Part 1: logit
------------------------------------------------------------------------------
Number of obs = 19386
LR chi2(2) = 2000.77
Prob > chi2 = 0.0000
Log likelihood = -8062.5899 Pseudo R2 = 0.1104
Part 2: regress_log
------------------------------------------------------------------------------
Number of obs = 15946
F( 2, 15943) = 1490.33
Prob > F = 0.0000
R-squared = 0.1575
Adj R-squared = 0.1574
Log likelihood = -29153.79 Root MSE = 1.5060
------------------------------------------------------------------------------
exp_tot | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
logit |
age | 0.047 0.001 33.81 0.000 0.045 0.050
female |
1 | 0.968 0.040 23.91 0.000 0.889 1.048
_cons | -0.871 0.060 -14.58 0.000 -0.988 -0.754
-------------+----------------------------------------------------------------
regress_log |
age | 0.036 0.001 52.82 0.000 0.034 0.037
female |
1 | 0.351 0.024 14.48 0.000 0.304 0.399
_cons | 5.329 0.037 142.80 0.000 5.256 5.402
------------------------------------------------------------------------------
与前文结果类似，年龄和性别的系数均在 1% 的水平上正显著，且 probit 模型与 logit 模型的 统计量也较为相似。并且，医疗支出的可能性与任何支出水平上的支出数量均随着年龄的增长的提高。与男性相比，女性更有可能花费至少 美元医疗支出，并且，在任何支出水平下，女性均比男性花费的更多。
随后，借助 margins 命令我们可以对医疗支出进行预测。
. margins, predict(duan) post
Predictive margins Number of obs = 19,386
Expression: twopm combined expected values, predict(duan)
------------------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | 4090.519 59.429 68.83 0.000 3974.041 4206.998
------------------------------------------------------------------------------
观察上表可知，预测的总医疗花费明显高于范例 3.1，数值约为每人每年 4091 美元。值得注意的是，当使用尺度再转换时，此处的 margins 命令无法提供正确的标准误。因此，我们必须借助非参数拔靴过程计算标准误与置信区间。
. program define Ey_boot, eclass
1. twopm exp_tot c.age i.female, firstpart(logit) secondpart(regress, log)
2. margins, predict(duan) nose post
3. end
. bootstrap _b, seed(14) reps(1000): Ey_boot
Predictive margins Number of obs = 19,386
Replications = 1,000
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | 4090.519 99.432 41.14 0.000 3895.635 4285.403
------------------------------------------------------------------------------
根据结果可知，拔靴过程计算所得的标准差几乎是增量法标准误的两倍大。因此，忽略尺度转换导致的不确定性将会严重低估标准差。与上述过程类似的是，借助 margins 命令计算边际效应时，仍需使用非参数拔靴过程进而获得正确的标准差。
. program define dydx_boot, eclass
1. twopm exp_tot c.age i.female, firstpart(logit) secondpart(regress, log)
2. margins, dydx(*) predict(duan) nose post
3. end
. bootstrap _b, seed(14) reps(1000): dydx_boot
Average marginal effects Number of obs = 19,386
Replications = 1,000
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age | 165.638 5.263 31.47 0.000 155.322 175.953
1.female | 1784.333 101.106 17.65 0.000 1586.170 1982.497
------------------------------------------------------------------------------
结果显示，年龄的边际效应约为 166 美元每年；同时，女性的边际效应约为 1784 美元。
3.3 总结
上述两种模型中不同的结果充分表明选择使用的模型确实很重要。但是，在没有进一步测试的情况下，尚不清楚哪个模型在统计意义上表现更好。正如 Duan（1983）所指出的，使用两部分模型可以产生实质性的改变，对 模型的重新转换方法也可以对结果产生很大的影响。上述原因均可能是我们示例中估算值之间差异的来源。
4. 结语
本篇推文中我们主要介绍了如何使用 twopm 这一命令估计两部分模型，进而为应对因变量为受限变量这一情况时提供了有力的武器。最后，值得一提的是， Tobit 模型与本文所介绍的两部分模型具有很大的相似性，但其也有很多不同，具体请参考周华林和李雪松 (2012)。
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh probit tobit heckman twopm, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata：广义Heckman两步法-gtsheckman Stata：两部模型实现命令-twopm Stata新命令：面板-LogitFE-ProbitFE 专题：回归分析 多层级Tobit模型及Stata应用 专题：面板数据 xtheckmanfe：面板Heckman模型的固定效应估计 专题：内生性-因果推断 Stata：样本选择偏误与两部模型-twopm-L121 Heckman 模型：你用对了吗？ Stata: 一文读懂 Tobit 模型 专题：交乘项-调节-中介 Logit-Probit中的交乘项及边际效应图示 专题：Probit-Logit Logit-Probit：非线性模型中交互项的边际效应解读 详解 Logit/Probit 模型中的 completely determined 问题 二元选择模型：Probit 还是 Logit？ Stata：二元Probit模型 动态 Probit 模型及 Stata 实现 专题：答疑-板书 FAQs答疑-2021寒假-Stata高级班-Day2-连玉君-面板门槛-Heckman-Tobit
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。