Stata负二项与零膨胀回归：处理过度离散与零值数据的完整指南

1. 项目概述：从泊松回归的局限说起

在实证研究的路上，尤其是处理计数数据时，泊松回归往往是我们的第一站。它的假设简洁明了：期望等于方差。但现实数据往往比教科书上的案例“调皮”得多。我处理过不少来自医学、社会学、经济学的数据集，比如一个社区全年的犯罪事件数、一家医院特定疾病的就诊人次、一个电商店铺的日投诉量。这些数据有一个共同点：它们都是非负整数，但方差常常远大于均值，我们称之为“过度离散”。当你用泊松回归拟合这类数据，得到的标准误会严重低估，导致你信心满满地认为发现了显著效应，实则可能只是模型误设带来的统计幻觉。

这时，负二项回归就该登场了。它通过引入一个额外的离散参数，优雅地放松了“均值=方差”的强假设，是处理过度离散计数数据的标准武器。但故事还没完。还有一种更棘手的情况：你的因变量里有一大堆零。比如，研究吸烟者每天的吸烟支数，很多人可能当天一支没抽（零值），而吸烟者则有一个正整数的分布。再比如，研究保险理赔次数，大部分保单持有人一年内可能零次理赔。当零值多到超出标准计数模型（泊松或负二项）的预测能力时，我们就遇到了“零膨胀”问题。此时，零膨胀模型，特别是零膨胀负二项回归，就成了解开数据谜团的关键钥匙。

今天，我们就深入Stata腹地，把nbreg和zinb这两个命令里里外外摸个透彻。这不仅仅是输入命令看结果，更是理解模型背后的逻辑、掌握Stata输出的每一行含义、并学会在复杂情境下（比如你搜索的“亚组分析”）灵活运用。我会结合多年实操中踩过的坑和总结的技巧，让你不仅能跑出回归，更能读懂数据在通过模型向你诉说的故事。

2. 核心模型原理与Stata命令逻辑拆解

2.1 负二项回归：泊松的“松绑”与离散参数alpha

为什么泊松回归会失灵？核心在于其方差与均值相等的强假设（Var(Y|X) = E(Y|X)）。负二项回归巧妙地引入了一个服从Gamma分布的误差项，使得条件方差成为条件均值的二次函数：Var(Y|X) = E(Y|X) + α*[E(Y|X)]^2。这里的α就是关键的超离散参数（overdispersion parameter）。当α = 0时，模型就退化成了泊松回归；α > 0则证实了过度离散的存在。

在Stata中，nbreg命令默认拟合的是NB2模型，即方差函数为上述的二次形式。这是最常见、最稳健的选择。命令基础语法很简单：

nbreg depvar [indepvars] [if] [in] [weight], options

但魔鬼藏在细节里。最重要的选项是dispersion()参数。虽然模型名为“负二项”，但Stata默认估计的是ln(α)，以保证其值非负。在输出结果中，你会看到一行/lnalpha的估计值及其标准误。真正的α需要通过di exp(_b[/lnalpha])来计算。许多新手会直接忽略这个/lnalpha，导致无法正确解读离散程度。

注意：nbreg默认使用最大似然估计。对于某些极端数据，α可能会非常大，提示可能存在零膨胀或其他结构性问题，这时就需要考虑zinb了。

2.2 零膨胀负二项回归：两个过程的混合

零膨胀模型的思想很直观：数据中的零来自两个不同的生成过程。

1. “必然零”过程：一个逻辑斯蒂（Logit）或概率（Probit）模型决定某个观测是否“必然”为零（例如，非吸烟者、无风险保单持有人）。这部分零无法用计数过程解释。
2. 计数过程：一个泊松或负二项模型，用于描述那些非“必然零”的观测（即可能取零也可能取正整数的观测）的计数分布。

因此，zinb命令实际上是在同时估计两个子模型：

• 膨胀模型：通常是Logit模型，预测“必然零”的概率。
• 计数模型：一个负二项回归模型，预测在非“必然零”状态下的计数期望。

其基本语法为：

zinb depvar [indepvars], inflate(varlist) [options]

这里的inflate()选项指定了哪些变量用于预测“必然零”的概率。这是模型设定中最需要理论思考的部分。你需要根据学科知识判断，哪些因素可能导致一个观测“根本不可能发生事件”。例如，在研究疾病发作次数时，inflate()里可以放入是否接种疫苗的变量，因为接种者可能“根本不可能”感染。

2.3 模型选择：如何决定用nbreg还是zinb？

这是一个实践性极强的问题，不能只看似然比检验。我通常遵循以下流程：

1. 初步诊断：先用poisson命令拟合，然后执行estat gof。如果卡方检验显著，表明泊松模型不合适，存在过度离散或零膨胀。
2. 检验过度离散：运行nbreg后，重点关注/lnalpha的估计值。对其进行假设检验（test _b[/lnalpha]=0），如果显著不为零，则支持负二项模型优于泊松模型。更直观的是看α的置信区间（通过nbreg, dispersion(mean)或事后计算）。
3. 检验零膨胀：这是关键。有两种常用方法：
   • Vuong检验：在zinb命令后使用vuong选项。该检验用于比较零膨胀模型与标准负二项模型。如果Vuong统计量显著为正，则支持零膨胀模型；显著为负则支持标准模型；不显著则两者难分优劣。
   • 计数拟合优度检验：使用countfit命令（需安装：ssc install countfit）。它可以同时比较泊松、负二项和零膨胀模型的拟合优度，提供非常直观的图表。
4. 理论依据：统计检验必须与理论结合。即使Vuong检验显著，你也必须能合理解释inflate()部分中变量的含义。如果找不到合理的变量来解释“必然零”过程，那么即使统计上显著，模型也可能缺乏实际意义。

3. 完整实操流程：从数据准备到结果解读

3.1 数据准备与探索性分析

在跑任何模型之前，彻底的描述性分析是必须的。假设我们有一个数据集health.dta，其中doc_visits表示一年内就诊次数，自变量有age,chronic（慢性病数量），insurance（是否有保险），gender等。

use health.dta, clear sum doc_visits tab doc_visits // 查看零值的比例 hist doc_visits, discrete freq // 绘制分布直方图

通过tab命令，你可能会发现doc_visits中零的比例高达40%。这是一个强烈的零膨胀信号。同时，计算方差与均值：

sum doc_visits, detail di r(Var)/r(mean)

如果比值远大于1（比如>1.5），则初步判断存在过度离散。

3.2 执行负二项回归

我们首先拟合一个标准的负二项模型。

nbreg doc_visits age chronic i.insurance gender, nolog

• nolog选项可以抑制迭代过程输出，让结果更清晰。
• i.insurance使用了因子变量语法，Stata会自动为分类变量生成虚拟变量。

结果解读要点：

1. 首先看模型整体的似然比检验（LR chi2）。它检验所有自变量系数是否联合为零。如果P值很小，说明模型整体显著。
2. 看各自变量的系数、标准误、Z值和P值。负二项回归的系数解释与泊松类似：exp(b)表示发生率比。例如，chronic的系数为0.3，则di exp(0.3) ≈ 1.35，意味着每增加一种慢性病，就诊次数的期望值将增加约35%。
3. 最关键的是看最底部的/lnalpha。运行test _b[/lnalpha]=0。如果拒绝原假设，则证实了过度离散的存在，使用负二项回归是合理的。记下alpha的值（di exp(_b[/lnalpha])），它量化了离散程度。

3.3 执行零膨胀负二项回归

基于理论，我们可能认为“没有保险”的人更可能因为费用问题而根本不去就诊（即“必然零”）。我们将insurance放入膨胀部分。

zinb doc_visits age chronic gender, inflate(insurance) vuong nolog

• inflate(insurance)指定了膨胀模型的自变量。
• vuong选项请求进行Vuong检验。

结果解读要点： Stata的输出分为上下两部分：

1. 上半部分：计数模型。解读方式与nbreg结果类似，系数表示在非“必然零”的群体中，自变量对就诊次数期望的影响。
2. 下半部分：膨胀模型（Logit）。这里的系数解释需要小心。系数为正，表示该变量增加“必然零”的概率。例如，insurance的系数若为正且显著，则表示有保险（假设insurance=1代表有保险）反而增加了成为“必然零”（即零次就诊）的对数发生比？这听起来不合常理。这里就体现了设定和编码的重要性。通常，我们可能认为insurance=0（无保险）才导致“必然零”。所以需要检查变量编码，或者系数应为负才符合直觉。exp(b)表示“必然零”的发生比。
3. Vuong检验：输出结果末尾会给出Vuong统计量。记住：显著为正支持ZINB，显著为负支持NB，不显著则无法判断。

3.4 边际效应与预测：让结果更直观

系数和发生比有时不够直观。margins命令可以计算平均边际效应或在特定值处的预测值，这对于向非专业受众解释结果至关重要。

预测期望计数：

* 计算所有观测在ZINB模型下的平均预测就诊次数 margins * 分别计算有保险和无保险群体的平均预测就诊次数 margins, over(insurance) * 绘制慢性病数量从0到5变化时，预测就诊次数的变化图（假设有保险） marginsplot, ytitle(“Predicted Doctor Visits”)

计算“必然零”的概率：

* 预测每个观测成为“必然零”的概率 predict pr_infl, pr sum pr_infl * 比较有保险和无保险群体的平均“必然零”概率 mean pr_infl, over(insurance)

这些预测值能让你更具体地理解模型含义，例如，“模型预测，没有保险的人群中，约有60%的人属于‘根本不会去就诊’的群体”。

4. 高级应用与疑难排解

4.1 如何进行亚组分析？

你搜索的“stata如何做亚组分析”是一个很实际的需求。对于nbreg或zinb，不建议简单地分样本回归然后比较系数，因为标准误可能不稳定，且难以进行正式的组间差异检验。更推荐的方法是使用交互项。

例如，想研究chronic对就诊次数的影响在gender间是否存在差异：

nbreg doc_visits age chronic##i.gender i.insurance, nolog

chronic##i.gender会自动生成chronic的主效应、gender的主效应以及它们的交互项。交互项的系数如果显著，就说明gender调节了chronic的影响。然后可以用margins来可视化这种调节效应：

margins gender, dydx(chronic) marginsplot, xdimension(gender)

这条margins命令会分别计算在男性和女性群体中，chronic增加一个单位对就诊次数的平均边际效应，并进行比较。

4.2 模型诊断与稳健性检验

1. 拟合优度：使用countfit命令（需安装）进行图形化比较。它会将实际数据的分布与模型预测的分布进行对比，一目了然。

   ssc install countfit quietly: zinb doc_visits age chronic gender, inflate(insurance) countfit doc_visits

2. 异常值检测：预测计数并与实际值比较，计算Pearson残差。

   predict yhat predict resid, pearson scatter resid yhat

   寻找残差绝对值过大的点，它们可能是模型拟合不佳的观测。
3. 稳健标准误：对于可能存在异方差或聚类结构的数据（如来自不同医院的患者），使用vce(robust)或vce(cluster cluster_var)选项来获得更稳健的标准误。

   nbreg doc_visits age chronic i.insurance, vce(cluster hospital_id)

4.3 常见报错与解决思路

• “initial values not feasible” 或 “convergence not achieved”：
  • 原因：模型过于复杂、初始值不佳、数据分离（特别是膨胀模型）。
  • 解决：
    1. 尝试from()选项提供初始值。可以先跑一个nbreg模型，然后用mat b = e(b)保存系数，在zinb中使用from(b)。
    2. 简化模型，特别是膨胀部分的变量。
    3. 检查膨胀部分的自变量是否存在完全预测零或非零的情况（数据分离）。
• Vuong检验结果为“NaN”或缺失：
  • 原因：通常发生在两个模型拟合结果非常接近，或某个模型拟合极差时。
  • 解决：优先依赖理论和其他拟合优度指标（如AIC/BIC）进行模型选择。estat ic命令可以输出信息准则。
• 系数符号与预期相反：
  • 原因：膨胀模型系数的解释是反直觉的。正系数意味着增加“必然零”的概率。务必厘清变量编码和业务逻辑。
  • 解决：使用margins命令直接计算关键变量对“必然零”概率的边际效应，这比解释系数更直接。

5. 结果呈现与报告撰写技巧

5.1 制作专业回归表格

手动整理结果效率低下且易错。推荐使用esttab命令（ssc install esttab）一键生成出版级表格。

* 分别估计泊松、负二项、零膨胀负二项模型 quietly: poisson doc_visits age chronic i.insurance gender estimates store Poisson quietly: nbreg doc_visits age chronic i.insurance gender estimates store NB quietly: zinb doc_visits age chronic gender, inflate(insurance) estimates store ZINB * 输出到屏幕，包含系数、标准误和显著性星号 esttab Poisson NB ZINB, b(%9.3f) se(%9.3f) star(* 0.1 ** 0.05 *** 0.01) /// stats(N ll alpha, fmt(%9.0f %9.1f %9.3f) labels(“N” “Log Likelihood” “Alpha”)) /// title(“Table 1: Comparison of Count Data Models”) * 输出到Excel文件 esttab Poisson NB ZINB using “results.xlsx”, replace /// b(%9.3f) se(%9.3f) star(* 0.1 ** 0.05 *** 0.01) /// stats(N ll alpha, fmt(%9.0f %9.1f %9.3f) labels(“N” “Log Likelihood” “Alpha”))

这张表格可以清晰展示不同模型下系数的变化、拟合优度（对数似然值）以及关键参数alpha，便于读者比较。

5.2 将暂元变量导出到文本文件

你搜索的“stata 将暂元变量导出到txt”是自动化报告和结果复现的好习惯。假设你想把关键的系数和标准误导出。

* 运行模型 zinb doc_visits age chronic gender, inflate(insurance) * 将关键结果存入暂元 local beta_age = _b[age] local se_age = _se[age] local p_age = 2*(1-normal(abs(_b[age]/_se[age]))) local alpha = exp(_b[/lnalpha]) * 打开一个文本文件并写入 file open myfile using “model_results.txt”, write replace file write myfile “ZINB Model Results” _n file write myfile “=================” _n _n file write myfile “Age Coefficient: `beta_age’ (SE: `se_age’, p: `p_age’)” _n file write myfile “Alpha (dispersion): `alpha’” _n file close myfile

这样，你就可以在后续的脚本或报告中自动调用这些结果。

5.3 可视化：让模型结果说话

除了表格，图形是更强大的沟通工具。

1. 预测概率图：展示不同chronic水平下，就诊次数为0, 1, 2, …的概率。

   quietly: zinb doc_visits age chronic gender, inflate(insurance) margins, at(chronic=(0(1)5)) predict(pr(0)) // 预测就诊0次的概率 marginsplot, title(“Probability of Zero Visits”) ytitle(“Probability”) recast(line)

2. 组间比较图：使用marginsplot绘制带有置信区间的边际效应图或预测值图，如前文亚组分析示例。

掌握nbreg和zinb，意味着你拥有了处理现实世界中复杂计数数据的两把利器。核心在于理解数据背后的故事：过多的零和过大的方差是数据在向你发出信号。通过系统的模型比较、严谨的诊断和深入的结果解读，你能让模型真正服务于研究问题，而非被复杂的输出表格所迷惑。每一次分析，从数据探索到模型诊断，再到清晰呈现，都是一个与数据对话的完整过程。