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 零膨胀负二项回归两个过程的混合零膨胀模型的思想很直观数据中的零来自两个不同的生成过程。“必然零”过程一个逻辑斯蒂Logit或概率Probit模型决定某个观测是否“必然”为零例如非吸烟者、无风险保单持有人。这部分零无法用计数过程解释。计数过程一个泊松或负二项模型用于描述那些非“必然零”的观测即可能取零也可能取正整数的观测的计数分布。因此zinb命令实际上是在同时估计两个子模型膨胀模型通常是Logit模型预测“必然零”的概率。计数模型一个负二项回归模型预测在非“必然零”状态下的计数期望。其基本语法为zinb depvar [indepvars], inflate(varlist) [options]这里的inflate()选项指定了哪些变量用于预测“必然零”的概率。这是模型设定中最需要理论思考的部分。你需要根据学科知识判断哪些因素可能导致一个观测“根本不可能发生事件”。例如在研究疾病发作次数时inflate()里可以放入是否接种疫苗的变量因为接种者可能“根本不可能”感染。2.3 模型选择如何决定用nbreg还是zinb这是一个实践性极强的问题不能只看似然比检验。我通常遵循以下流程初步诊断先用poisson命令拟合然后执行estat gof。如果卡方检验显著表明泊松模型不合适存在过度离散或零膨胀。检验过度离散运行nbreg后重点关注/lnalpha的估计值。对其进行假设检验test _b[/lnalpha]0如果显著不为零则支持负二项模型优于泊松模型。更直观的是看α的置信区间通过nbreg, dispersion(mean)或事后计算。检验零膨胀这是关键。有两种常用方法Vuong检验在zinb命令后使用vuong选项。该检验用于比较零膨胀模型与标准负二项模型。如果Vuong统计量显著为正则支持零膨胀模型显著为负则支持标准模型不显著则两者难分优劣。计数拟合优度检验使用countfit命令需安装ssc install countfit。它可以同时比较泊松、负二项和零膨胀模型的拟合优度提供非常直观的图表。理论依据统计检验必须与理论结合。即使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, nolognolog选项可以抑制迭代过程输出让结果更清晰。i.insurance使用了因子变量语法Stata会自动为分类变量生成虚拟变量。结果解读要点首先看模型整体的似然比检验LR chi2。它检验所有自变量系数是否联合为零。如果P值很小说明模型整体显著。看各自变量的系数、标准误、Z值和P值。负二项回归的系数解释与泊松类似exp(b)表示发生率比。例如chronic的系数为0.3则di exp(0.3) ≈ 1.35意味着每增加一种慢性病就诊次数的期望值将增加约35%。最关键的是看最底部的/lnalpha。运行test _b[/lnalpha]0。如果拒绝原假设则证实了过度离散的存在使用负二项回归是合理的。记下alpha的值di exp(_b[/lnalpha])它量化了离散程度。3.3 执行零膨胀负二项回归基于理论我们可能认为“没有保险”的人更可能因为费用问题而根本不去就诊即“必然零”。我们将insurance放入膨胀部分。zinb doc_visits age chronic gender, inflate(insurance) vuong nologinflate(insurance)指定了膨胀模型的自变量。vuong选项请求进行Vuong检验。结果解读要点 Stata的输出分为上下两部分上半部分计数模型。解读方式与nbreg结果类似系数表示在非“必然零”的群体中自变量对就诊次数期望的影响。下半部分膨胀模型Logit。这里的系数解释需要小心。系数为正表示该变量增加“必然零”的概率。例如insurance的系数若为正且显著则表示有保险假设insurance1代表有保险反而增加了成为“必然零”即零次就诊的对数发生比这听起来不合常理。这里就体现了设定和编码的重要性。通常我们可能认为insurance0无保险才导致“必然零”。所以需要检查变量编码或者系数应为负才符合直觉。exp(b)表示“必然零”的发生比。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, nologchronic##i.gender会自动生成chronic的主效应、gender的主效应以及它们的交互项。交互项的系数如果显著就说明gender调节了chronic的影响。然后可以用margins来可视化这种调节效应margins gender, dydx(chronic) marginsplot, xdimension(gender)这条margins命令会分别计算在男性和女性群体中chronic增加一个单位对就诊次数的平均边际效应并进行比较。4.2 模型诊断与稳健性检验拟合优度使用countfit命令需安装进行图形化比较。它会将实际数据的分布与模型预测的分布进行对比一目了然。ssc install countfit quietly: zinb doc_visits age chronic gender, inflate(insurance) countfit doc_visits异常值检测预测计数并与实际值比较计算Pearson残差。predict yhat predict resid, pearson scatter resid yhat寻找残差绝对值过大的点它们可能是模型拟合不佳的观测。稳健标准误对于可能存在异方差或聚类结构的数据如来自不同医院的患者使用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”原因模型过于复杂、初始值不佳、数据分离特别是膨胀模型。解决尝试from()选项提供初始值。可以先跑一个nbreg模型然后用mat b e(b)保存系数在zinb中使用from(b)。简化模型特别是膨胀部分的变量。检查膨胀部分的自变量是否存在完全预测零或非零的情况数据分离。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 可视化让模型结果说话除了表格图形是更强大的沟通工具。预测概率图展示不同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)组间比较图使用marginsplot绘制带有置信区间的边际效应图或预测值图如前文亚组分析示例。掌握nbreg和zinb意味着你拥有了处理现实世界中复杂计数数据的两把利器。核心在于理解数据背后的故事过多的零和过大的方差是数据在向你发出信号。通过系统的模型比较、严谨的诊断和深入的结果解读你能让模型真正服务于研究问题而非被复杂的输出表格所迷惑。每一次分析从数据探索到模型诊断再到清晰呈现都是一个与数据对话的完整过程。
