Stata：一文读懂两部模型-twopm

发布时间：2023-03-08 阅读 384

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

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

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

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：王恩泽 (武汉大学经济与管理学院)
邮箱：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) $y_{i}$ ≥ 0，(2) $y_{i}$ = 0。换句话说，该类变量往往由大量的 $0$ 值与大于 $0$ 的值组成。由于大量 $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 对 $l n (y)$ 进行回归，该程序能够提供以 $y$ 尺度 ( 原始尺度 ) 为基础预测的估计值。
第五，该命令会在考虑模型的两个部分、复杂的调查设计、及基于增量方法 (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：为模型的每个部分创建得分变量。由于模型第二部分的得分只有在估计 $y > 0$ 分样本时有意义，故计算第二部分的得分变量时样本被自动限于估计的分样本。
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 模型，在第二步中对对数形式的 $y$ 使用了 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 模型的 $z$ 统计量也较为相似。并且，医疗支出的可能性与任何支出水平上的支出数量均随着年龄的增长的提高。与男性相比，女性更有可能花费至少 $1$ 美元医疗支出，并且，在任何支出水平下，女性均比男性花费的更多。

随后，借助 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）所指出的，使用两部分模型可以产生实质性的改变，对 $l n (y)$ 模型的重新转换方法也可以对结果产生很大的影响。上述原因均可能是我们示例中估算值之间差异的来源。

4. 结语

本篇推文中我们主要介绍了如何使用 twopm 这一命令估计两部分模型，进而为应对因变量为受限变量这一情况时提供了有力的武器。最后，值得一提的是， Tobit 模型与本文所介绍的两部分模型具有很大的相似性，但其也有很多不同，具体请参考周华林和李雪松 (2012)。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh probit tobit heckman twopm, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]