专栏名称: 连享会

连玉君老师团队分享，主页：lianxh.cn。白话计量，代码实操；学术路上，与君同行。

偏态分布数据的回归模型选择

连享会 · 公众号 · · 2024-06-22 22:52

正文

👇 连享会 · 推文导航 | www.lianxh.cn

🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令
📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归
⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析
🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
🔨 工具：工具软件 | Markdown | Python-R-Stata
🎧 课程：公开课-直播 | 计量专题 | 关于连享会

连享会课程 · 2024 政策分析专题

作者：刘潍嘉 (兰州大学)
邮箱：[email protected]

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

1. 成本数据的特点

成本数据的分布通常呈现偏态，伴随着长右尾特点，同时也包含相当比例的零值。尽管普通最小二乘法 (OLS) 可以很好地预测平均成本，但它并不能很好地拟合这种分布特性，因此在这种情况下不能直接应用 OLS。本文将以医疗成本为例，探讨偏态分布数据作为被解释变量时的回归模型选择与应用。

本文采用 2017 年 Medical Expenditure Panel Survey (MEPS) 的数据进行分析，旨在评估被诊断为高血压的家庭受访者的平均总医疗支出。通过比较不同回归模型的结果，本文探讨了模型选择的问题。在分析过程中，我们控制了多个样本特征，包括年龄、性别、种族、民族、贫困状况、婚姻状况以及人口普查地区。

本文的前三章使用了 Stata 进行数据导入和分析，而第四章则应用 R 进行案例分析。全文代码整理在本文附录中。

. * 导入数据
. clear all
. lxhget limited_data.csv, replace 
. import delimited using limited_data.csv, clear

. * 分析中仅包括18岁及以上的受访者
. keep if age17x >= 18

. * 清洗数据
. gen female = .
. replace female = 0 if sex == "1 MALE"
. replace female = 1 if sex == "2 FEMALE"
. label define female_lbl 0 "Male" 1 "Fenale"
. label values female female_lbl
. tab female, m
. gen race = .
. replace race = 0 if racev2x == "1 WHITE - NO OTHER RACE REPORTED"
. replace race = 1 if racev2x == "2 BLACK - NO OTHER RACE REPORTED"
. replace race = 2 if racev2x == "3 AMER INDIAN/ALASKA NATIVE-NO OTHER RACE"
. replace race = 3 if racev2x == "4 ASIAN INDIAN - NO OTHER RACE REPORTED"
. replace race = 3 if racev2x == "5 CHINESE - NO OTHER RACE REPORTED"
. replace race = 3 if racev2x == "6 FILIPINO - NO OTHER RACE REPORTED"
. replace race = 3 if racev2x == "10 OTH ASIAN/NATV HAWAIIAN/PACFC ISL-NO OTH"
. replace race = 4 if racev2x == "12 MULTIPLE RACES REPORTED"
. label define race_lbl 0 "White" 1 "Black" 2 "American Indian/Alaska Native" ///
> 3 "Asian" 4 "Multiple races reported"
. label values race race_lbl
. tab race, m
. gen hispanic = .
. replace hispanic = 0 if hispanx == "2 NOT HISPANIC"
. replace hispanic = 1 if hispanx == "1 HISPANIC"
. tab hispanic, m
. gen marital_status = .
. replace marital_status = 0 if marry17x == "5 NEVER MARRIED"
. replace marital_status = 1 if marry17x == "1 MARRIED"
. replace marital_status = 2 if marry17x == "2 WIDOWED"
. replace marital_status = 3 if marry17x == "3 DIVORCED"
. replace marital_status = 4 if marry17x == "4 SEPARATED"
. replace marital_status = 5 if marry17x == "-7 REFUSED"
. replace marital_status = 5 if marry17x == "-8 DK"
. label define marital_lbl 0 "Never married" 1 "Married" 2 "Widowed" 3 "Divorced" ///
> 4 "Separated" 5 "Unknown/Refused"
. label values marital_status marital_lbl
. tab marital_status, m
. gen poverty_status = .
. replace poverty_status = 0 if povcat17 == "1 POOR/NEGATIVE"
. replace poverty_status = 1 if povcat17 == "2 NEAR POOR"
. replace poverty_status = 2 if povcat17 == "3 LOW INCOME"
. replace poverty_status = 3 if povcat17 == "4 MIDDLE INCOME"
. replace poverty_status = 4 if povcat17 == "5 HIGH INCOME"
. label define poverty_lbl 0 "Poor" 1 "Near poor" 2 "Low income" ///
> 3 "Middle income" 4 "High income"
. label values poverty_status poverty_lbl
. tab poverty_status, m 
. gen region = .
. replace region = 0 if region17 == "1 NORTHEAST"
. replace region = 1 if region17 == "2 MIDWEST"
. replace region = 2 if region17 == "3 SOUTH"
. replace region = 3 if region17 == "4 WEST"
. label define region_lbl 0 "Northeast" 1 "Midwest" 2 "South" 3 "West"
. label values region region_lbl
. tab region, m
        
. * 重命名变量
. drop region17 sex racev2x hispanx marry17x povcat17
. rename region region17
. rename female sex
. rename race racev2x
. rename hispanic hispanx
. rename marital_status marry17x
. rename poverty_status povcat17
        
. * 可视化以初步观察数据分布情况
. histogram totexp17, bin(50) percent normal
. kdensity totexp17

. * 偏度分析
. summarize totexp17, detail

                          totexp17
-------------------------------------------------------------
      Percentiles      Smallest
 1%            0              0
 5%            0              0
10%          161              0       Obs               7,872
25%        953.5              0       Sum of wgt.       7,872

50%         3517                      Mean            10625.1
                        Largest       Std. dev.       23462.3
75%      10549.5         474178
90%        26858         499286       Variance       5.50e+08
95%        43280         506064       Skewness       8.581631
99%       105557         552898       Kurtosis       136.4949

. sktest totexp17 /* 偏度检验 */

上图分别展示了医疗支出（成本）的直方图和核密度图，可以观察到数据分布的两个明显特点：

在接近 0 的区域存在大量观测值。
数据分布呈现长右尾现象，有一小部分观测值明显偏离主要数据分布，且数值相对较大。

2. 各模型的对比分析

在对比分析中，我们将分别应用以下模型 (OLS，Log-OLS，Log-OLS with smearing，GLM，two-part models) 进行分析，然后对比各模型的拟合优劣，以及拟合结果的统计特征。在本文中，模型优劣（Goodness of Fit tests，简称 GOF）主要通过以下三个检验进行对比分析：

Pearson 相关性检验：检验成本拟合值和残差之间的相关性。在此检验中，我们期望拟合值和残差之间不存在任何相关性。
Pregibon’s Link 检验：检验模型是否被正确设定。在此检验中，我们对拟合值进行平方，并用拟合值和其平方作为解释变量再次进行回归。如果拟合值的平方与成本没有关联，则模型设定被认定为正确。
The modified Hosmer-Lemeshow 检验：检验平均残差是否等于零。通常使用 10 个十分位数的拟合值进行此检验，如果检验结果在统计上不显著，则模型设定被认定为正确。

2.1 OLS

. ** 模型 1: OLS
. **** 模型 1: OLS
. reg totexp17 age17x sex racev2x hispanx marry17x povcat17 region17 /* OLS model */ 
. predict yhat /* 得到拟合值 */
. predict error, resid /* 得到残差 */
. graph twoway scatter error yhat /* 绘制散点图 */
. rvfplot, yline(0) /* 绘制散点图同样也可以直接应用这个命令 */
. summarize yhat, detail
. summarize totexp17 yhat
. histogram yhat, bin(50) percent normal
. estat szr age17x sex racev2x hispanx marry17x povcat17 region17  /* Szroeter's 异方差性检验 */

. **** GOF tests
. qui reg totexp17 age17x sex racev2x hispanx marry17x povcat17 region17 
. predict xb, xb
. gen mu_ols = xb
. gen res_ols = totexp17 - mu_ols
. * Pearson corr
. pwcorr res_ols mu_ols, sig
. * Link test
. gen xb2 = xb^2
. reg totexp17 xb xb2
. qui reg totexp17 age17x sex racev2x hispanx marry17x povcat17 region17 
. linktest
. * H-L test
. xtile xbtile = xb, nq(10)
. qui tab xbtile, gen(xbt)
. reg res_ols xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10, nocons
. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10
. drop xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10 xbtile xb xb2

首先对 OLS 模型进行分析，其中解释变量包括年龄、性别、种族、民族、贫困状况、婚姻状况以及人口普查地区。

线性模型提供了对系数的简单解释。
然而，由于数据分布有较高的偏度，尾部的差异会对回归结果产生很大影响。
由于 Y 的非线性，进而产生有偏估计。
并且可能因异方差性产生无效的标准误。

我们应用 OLS 进行回归后，绘制了残差与拟合值的散点图。可以看出，残差分布随着拟合值的增大变得更加分散，这可能表明存在异方差问题。

然后，我们对成本的拟合值和原始值的数据特征进行了讨论。结果显示，虽然拟合值与原始值的平均值相同，但它们的标准差、最小值和最大值各不相同。尤其是拟合值中出现了负值，对于医疗成本来说，这显然是不合理的。

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
        yhat |      7,872     10625.1    3213.576  -462.9974   19580.85

我们进一步绘制了拟合值的直方图，发现拟合值的分布发生了变化，呈现出正态分布的特征。

然后，我们对 OLS 回归结果进行了之前提到的三项检验，得到了以下的 GOF 检验结果：

Pearson 相关系数：P = 1.000
Pregibon 的联合检验：P = 0.004
经过修正的 Hosmer-Lemeshow 检验：P = 0.539

检验 1 表明残差与拟合值之间没有相关性，这是我们期望的结果。然而，检验 2 显示拟合值的平方与结果之间存在显著相关性，这是我们不期望看到的。最后，检验 3 表明平均残差与 0 无显著差异，这也是我们期望的结果。因此，尽管检验 1 和检验 3 符合我们的期望，但检验 2 的结果仍然表明我们的模型设定可能存在问题。

2.2 Log - OLS

. ** 模型 2: Log-OLS
. gen lntptexp = ln(totexp17)
. histogram lntptexp, bin(50) percent normal
. summarize lntptexp, detail

. **** 模型 2: Log-OLS
. reg lntptexp age17x sex racev2x hispanx marry17x povcat17 region17
. predict lh_yhat, xb
. gen exp_lnyhat = exp(lh_yhat)
. summarize exp_lnyhat, detail
. summarize totexp17 exp_lnyhat
 
. * Pearson corr
. gen res_lny = totexp17 - exp_lnyhat
. pwcorr res_lny exp_lnyhat, sig
 
. * Link test
. gen xb2 = lh_yhat^2
. reg lntptexp lh_yhat xb2, robust

. * H-L test
. xtile xbtile = lh_yhat, nq(10)
. qui tab xbtile, gen(xbt)
. reg res_lny xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10, nocons robust
. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10
. drop xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10 xbtile xb2

对成本数据进行对数转换可以降低偏度 (如下图所示)，因此我们在进行对数转换后再次应用 OLS 进行分析。

然而，必须指出，对数转换存在一定的局限性，并且其回归系数的含义也会随之改变，具体内容请参见如下推文。

秦范, 2021, 取对数！取对数？, 连享会 No.615.
毕英睿, 2023, 取对数：如何应对零值和负数, 连享会 No.1234.
吴梦萱, 2024, 纠结！DID 中取对数还是不取对数？论文推介, 连享会 No.1340.

Log - OLS 模型的成本原始值与拟合值的描述性统计如下，另外，对数化后的数据分布情况如上图所示，基本呈正态分布。

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
  exp_lnyhat |      7,872    3914.653    1966.329   483.2885   12193.77

遵循相同的方法，我们对 Log-OLS 模型的结果进行了三项检验，GOF 检验的结果如下：

Pearson 相关性检验：P = 0.000
Pregibon 的链路检验：P = 0.279
修改后的 Hosmer-Lemeshow 检验：P = 0.000

检验 1 表明残差与拟合值之间存在显著相关性，这意味着可能存在异方差性，而这是我们不希望看到的结果。检验 2 显示拟合值的平方与结果之间没有显著相关性，但是检验 3 显示平均残差与 0 之间存在显著差异。因此，我们的模型设定仍存在问题。

2.3 Log - OLS smearing estimator


. ** 模型 3: Log-OLS w/smearing
. *** 模型 3: Log-OLS w/smearing
. reg lntptexp age17x sex racev2x hispanx marry17x povcat17
. predict xb, xb

. * Smearing estimator
. gen smr = exp(lntptexp - lh_yhat)
. summarize smr
. gen smear = r(mean) 
. gen mu = exp(lh_yhat) * smear
. gen mu_lols = exp(lh_yhat) * smear
. gen res_lols = totexp17 - mu_lols
. summarize totexp17 mu 
. summarize mu, detail

. * Park test - 检验异方差
. gen res = (lntptexp - xb)^2
. glm res age17x sex racev2x hispanx marry17x povcat17, family(gamma) link(log) robust
. test

. * Pearson corr
. pwcorr res_lols mu, sig

. * Link test
. gen xb2 = xb^2
. reg lntptexp xb xb2, robust

. * H-L test
. xtile xbtile = xb, nq(10)
. qui tab xbtile, gen(xbt)
. reg res_lols xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10, nocons robust
. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10
. drop xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10 xbtile xb2

在章节 2.2 中，我们应用 Log-OLS 模型对医疗成本进行了估计。然而，如数据的描述性统计所示，我们得到的是对数化后的拟合值。那么，我们如何得到实际的成本呢？是否可以通过对拟合值取反对数来直接转换？

对于模型，我们尝试通过取反对数来获得拟合值。

只有当我们假设时，我们才能得到。但是，当时，能否认为呢？即是否成立？通过以下一个简单实例可以得出，答案是否定的。残差的反对数的期望值不等于残差期望值的反对数。

当并且时：

因此，。

Duan (1983) 提出的 smearing estimator 有效地解决了这一问题，通过以下公式进行转换。

经 smearing estimator 调整后的回归，转换后的拟合值与原始值的统计描述如下：

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
          mu |      7,872    12588.82    6323.361   1554.169   39212.99

对 Log-OLS smearing estimator 结果进行了三项检验，GOF 检验结果如下：

Pearson 相关性检验：P = 0.000
Pregibon’s Link 检验：P = 0.267
修改后的 Hosmer-Lemeshow 检验：P = 0.000

三项检验结果与 Log-OLS 模型的检验结果相似，因此模型设定仍然存在问题。

2.4 Generalized Linear Model (GLM)

. ** MODEL 4: GLM (gamma)
. glm totexp17 age17x sex racev2x hispanx marry17x povcat17, family(gamma) link(log)
. predict glm_1
. summarize glm_1, detail
. histogram glm_1, bin(50) percent normal
. drop xb
. predict xb, xb
. gen mu_glm = exp(xb)
. gen res_glm = totexp17 - mu_glm
. summarize totexp17 mu_glm 

. * Pearson corr
. pwcorr res_glm mu_glm, sig

. * Link test
. gen xb2 = xb^2
. reg totexp17 xb xb2, robust

. * H-L test
. xtile xbtile = xb, nq(10)
. qui tab xbtile, gen(xbt)
. reg res_glm xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10, nocons robust
. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10
. drop xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10 xbtile xb xb2

广义线性模型 (GLM) 是对 OLS 的扩展。GLM 假设我们关注的随机变量的分布可以通过链接函数 (link function) 与系统性效应 (即非随机效应) 建立关联，从而解释其相关性。它允许每个观测值的方差大小与其拟合值相关，从而扩展了线性回归的应用范围。

在此，GLM 使用链接函数。我们变换的不是原始的，而是。在模型设定中，Family 的选择基于和之间的关系，本文选择了 Gamma 分布。而链接函数的选择则基于 Pregibon’s link test。经 GLS 回归得到后的拟合值数据分布如下：

通过 GLS 回归得到的拟合值与原始值的统计描述如下：

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
      mu_glm |      7,872    10639.42    3261.605   3534.286   21606.15

对于 GLS 结果进行三项检验，GOF 检验结果如下。

Pearson correlation：P = 0.426
Pregibon’s Link test：P = 0.534
The modified Hosmer-Lemeshow test：P = 0.910

检验 1 表明残差与拟合值之间没有显著的相关性。检验 2 显示拟合值的平方与结果之间也没有显著的相关性。检验 3 表明平均残差与 0 之间没有显著差异。因此，经过 GLS 回归和上述检验结果表明模型没有设定错误，是一个相对可行的方案。

2.5 Two-Part model

. ** MODEL 5: Two-part model
. twopm totexp17 age17x sex racev2x hispanx marry17x povcat17, firstpart(logit) ///
> secondpart(glm, family(gamma) link(log))
. predict twopm_xb
. summarize twopm_xb, detail
. histogram twopm_xb, bin(50) percent normal
. gen res_twopm = totexp17 - twopm_xb
. summarize totexp17 mu mu_glm twopm_xb

. * Pearson corr
. pwcorr res_twopm twopm_xb, sig

. * Link test
. gen xb2 = twopm_xb^2
. reg totexp17 twopm_xb xb2, robust
. xtile xbtile = twopm_xb, nq(10)
. qui tab xbtile, gen(xbt)
. reg res_twopm xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10, nocons robust
. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10
. drop xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9 xbt10 xbtile

我们关注的医疗成本数据中包含大量的 0 值。鉴于这些大量的 0 值，我们需要使用两部分模型 (Two-part models) 来处理这类数据。两部分模型考虑到了大量零值的存在以及成本分布的非参数属性。

两部分模型的第一部分使用 Logit、Probit 等方法来预测事件是否发生，即被解释变量是否为 0。这一部分的目标是评估被解释变量是否大于零，即发生医疗成本的概率。在第二部分中，对于那些在第一部分中被判断为发生医疗成本的个体，使用 GLM 来对这些样本进行拟合。因此，的期望值考虑了样本成本是否为非零的情况。

通过 Two-Part 模型回归得到的拟合值数据分布如下：

经 Two-Part model 回归得到后的拟合值与原始值统计性描述如下：

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
    twopm_xb |      7,872    10625.31    3160.112   2796.703   20732.22

对于 Two-Part model 结果进行三项检验，GOF 检验结果如下。

Pearson correlation：P = 0.782
Pregibon’s Link test：P = 0.129
The modified Hosmer-Lemeshow test：P = 0.788

检验 1 表明残差与拟合值之间没有显著的相关性。检验 2 显示拟合值的平方与结果之间不存在显著的相关性。检验 3 显示平均残差与 0 没有显著差异。因此，经过 Two-Part 模型回归和上述检验结果表明，该模型没有设定错误，是相对可行的方案。

3. 结论

最后，我们列出了成本的原始值与五个模型拟合值的描述性统计并进行对比，发现应用 OLS 估计 (yhat) 得到的结果出现了负值，这显然与现实不符。而 Log-OLS 结果 (exp_lnyhat) 的再转换会引入偏差，使得难以分析其现实意义。对于 Log-OLS smearing estimator、GLM 以及 Two-Part model，三者的估计结果均能够令人接受，尤其是后两项估计，其拟合值的均值与原始值的均值最为接近。

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    totexp17 |      7,872     10625.1     23462.3          0     552898
        yhat |      7,872     10625.1    3213.576  -462.9974   19580.85
  exp_lnyhat |      7,872    3914.653    1966.329   483.2885   12193.77
          mu |      7,872    12588.82    6323.361   1554.169   39212.99
      mu_glm |      7,872    10639.42    3261.605   3534.286   21606.15
    twopm_xb |      7,872    10625.31    3160.112   2796.703   20732.22

进一步地，我们将三项检验结果汇总如下，GLM 和 Two-Part 模型均通过了所有三项检验。因此，我们得出结论，对于类似本文偏态数据分布的情况，应优先考虑应用 GLM 和 Two-Part 模型。

	OLS	Log - OLS	Log - OLS smearing	GLM	Two-Part model
Pearson correlation	1.000	0.000	0.000	0.426	0.782
Pregibon’s Link test	0.004	0.279	0.267	0.534	0.129
The modified Hosmer-Lemeshow test	0.539	0.000	0.000	0.910	0.788

4. 两部模型在 R 中的应用案例

为了探讨如何使用两部模型，我们将利用 2017 年医疗支出面板调查 (MEPS) 数据中高血压受访者的总医疗支出。具体而言，我们希望估算 2017 年美国高血压患者的平均医疗支出。我们将重点关注美国 (US) 患有高血压的成年人 (年龄 >= 18 岁)。

> # 导入数据
> data1 "https://file.lianxh.cn/data/l/limited_data.csv")
> # 限制数据为 >= 18岁的成年样本
> data1 = "18"), ]
> hist(data1$totexp17)

> summary(data1$totexp17)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     0.0    953.8   3517.0  10625.1  10547.8 552898.0

观察数据分布和描述性统计，发现：

平均医疗保健支出为每人 $10,625。
有很多人医疗保健支出为零或接近零。
医疗保健支出分布的长右尾，这是因为存在少数医疗保健支出较大的样本。

要构建两部模型，我们需要将 twopartm 包安装并加载到 R 中。

观察数据分布和描述性统计结果后，我们可以得出以下结论：

平均医疗保健支出为每人 $10,625。
很多人的医疗保健支出接近零或为零。
医疗保健支出呈现长尾分布，即存在一小部分个体的医疗保健支出较高。

为了构建两部模型，我们需要先安装并加载 twopartm 包到 R 环境中。

# 安装twopartm包
install.packages("twopartm")
# 加载twopartm包
> library("twopartm")

我们使用的主要函数是 tpm。通常情况下，两部模型的第一部分和第二部分会使用相同的变量，但有时我们也会希望为这两部分设定不同的模型。在第一部分的设定中，我们主要使用 formula_part1，而第二部分的设定则使用formula_part2。在本案例中，我们将在第一部分中使用 logit 模型。第二部分将采用 gamma 模型，以捕获从零到无穷大的成本分布范围（非负值）。

> # Two-part model
> tpm.model1 +                   data = data1, link_part1 = "logit", family_part2 = Gamma(link = "log"))
> summary(tpm.model1)

$Firstpart.model
Call:
glm(formula = nonzero ~ age17x + sex + racev2x + hispanx + marry17x + 
    povcat17, family = binomial(link = "logit"), data = data1)


$Secondpart.model
Call:
glm(formula = totexp17 ~ age17x + sex + racev2x + hispanx + marry17x + 
    povcat17, family = Gamma(link = "log"), data = data1)

在输出的结果中，我们获得了两部模型第一部分和第二部分的详细信息。然而，这些结果并未直接提供 2017 年美国高血压患者的平均医疗总支出。为了获取这一关键信息，我们需运用 predict 函数来估计样本的平均医疗保健总支出。此外，我们还需计算其标准误差。通过在 predict 函数中设置 se.fit = TRUE 选项，我们能够运用 Delta 方法来估计这一标准误差。

> # 生成每个样本的总成本估计值
> predict1 TRUE)
> summary(predict1$fit)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2797    8365   10419   10625   12600   20732 

> summary(predict1$se.fit)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  345.0   542.7   692.3   770.7  936.7  2533.9

根据 predict 函数的结果，2017 年美国成年高血压患者的平均医疗支出为 10,625 美元，标准误差为 770.7 美元。我们可以使用以下公式计算 95% 置信区间 (CI)，从而得到下限 (LL) 和上限 (UL) 来表示 95% CI 的范围：

> predict1$ci 1.96*predict1$se.fit
> predict1$ll > predict1$ul > summary(predict1$ll)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1834    7142    9022    9115   10921   17320

> summary(predict1$ul)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3473    9575   11740   12136   14330   24479

根据计算结果，95% 置信区间介于 9,115 美元到 12,136 美元之间。得到这一结果后，我们可以进行拟合优度 (GOF) 测试，以评估我们的模型设定是否正确。我们将使用前文提到的三种 GOF 测试来对模型进行检验。

> # GOF tests
> res ## 估计残差
> summary(res)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-20492.7  -8725.3  -5974.2     -0.2   -468.3 542082.9

> # Pearson correlation
> pwcorr > pwcorr

  Pearson's product-moment correlation
 data:  predict1$fit and res
 t = -0.27607, df = 7870, p-value = 0.7825
 alternative hypothesis: true correlation is not equal to 0
 95 percent confidence interval:
  -0.02520129  0.01898051
 sample estimates:
          cor 
 -0.003111908

> # Pregibon's link test
> xb2 ## 将预测值平方
> linear.model +                     family = gaussian(link = "identity"))
> summary(linear.model)

Call:
glm(formula = data1$totexp17 ~ predict1$fit + xb2, family = gaussian(link = "identity"))
Deviance Residuals: 
    Min      1Q  Median      3Q     Max  
 -17864   -8937   -5961    -479  541759  
 Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
 (Intercept)  -3.522e+03  2.432e+03  -1.448 0.147620    
 predict1$fit  1.704e+00  4.424e-01   3.851 0.000118 ***
 xb2          -3.219e-05  1.925e-05  -1.672 0.094512 .  
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 (Dispersion parameter for gaussian family taken to be 540891244)
     Null deviance: 4.3328e+12  on 7871  degrees of freedom
 Residual deviance: 4.2563e+12  on 7869  degrees of freedom
 AIC: 180641
 Number of Fisher Scoring iterations: 2

> # Modified Hosmer-Lemeshow test with 10 deciles
> library("ResourceSelection")  ## 加载包
> hoslem.test(data1$totexp17, predict1$fit, g = 10)

  Hosmer and Lemeshow goodness of fit (GOF) test
 
 data:  data1$totexp17, predict1$fit
 X-squared = -30.975, df = 8, p-value = 1

根据拟合优度 (GOF) 测试的结果，验证显示我们的两部分模型设定基本正确，从而确认其对 2017 年美国高血压患者总医疗支出的估计是可靠的。

5. 参考资料

Mark Bounthavong. (2023). Cost as a dependent variable -Link-
Paul G. Barnett. (2017). Econometrics Course: Cost as the Dependent Variable (I) -Link-
Paul G. Barnett. (2017). Econometrics Course: Cost as the Dependent Variable (II) -Link-
Duan, N. (1983). Smearing estimate: A nonparametric retransformation method. Journal of the American Statistical Association, 78(383), 605–610. -Link-
Cost as a dependent variable: Application of the two-part model -Link-
Github - Cost as a dependent variable -Link-

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 正态对数两部, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题：倍分法DID

吴梦萱, 2024, 纠结！DID 中取对数还是不取对数？论文推介, 连享会 No.1340.

专题：Stata命令

张少鹏, 2021, Stata：两部模型实现命令-twopm, 连享会 No.742.

专题：内生性-因果推断

张少鹏, 2021, Stata：样本选择偏误与两部模型-twopm-L121, 连享会 No.714.

专题：数据处理

毕英睿, 2023, 取对数：如何应对零值和负数, 连享会 No.1234.
王恩泽, 2023, Stata：一文读懂两部模型-twopm, 连享会 No.1170.

专题：Stata教程

王硕, 2023, Stata：正态转换的五种方法, 连享会 No.1237.

专题：回归分析

秦范, 2021, 取对数！取对数？, 连享会 No.615.
连玉君, 杨柳, 2020, Stata: 边际效应分析, 连享会 No.64.
陈卓然, 2022, Stata：多元正态t截断分布的命令, 连享会 No.883.
陈宇灿, 2022, Stata：线性与对数线性函数形式选择-imfreg, 连享会 No.1152.

专题：Probit-Logit

陈宇灿, 2023, Stata：线性与对数线性函数形式选择-imfreg, 连享会 No.1153.