很多人第一次接触 Stata是因为导师或同事丢过来一个 do 文件说照着跑一遍就行。我自己当年也是这么入门的一份几十个变量的问卷数据一张表格要填十几个模型的系数用表格软件硬算根本算不出来。后来项目做得多了从横截面问卷到面板数据从只盯 p 值到认真处理内生性和异质性多元回归分析在Stata里几乎成了我每周都要碰的东西。这篇东西不打算复述教科书只讲实操里踩过、绕过、最后固定下来的一整套流程装完之后先做什么、regress输出那张表每一列到底在说什么、稳健标准误什么时候必须加、亚组分析该分组回归还是做交互、报错信息怎么读以及一堆看起来不起眼但能省半天时间的细节命令。刚上手的学生能照着走完一遍做了几年想回头补基础的从业者也能从里面挑到能直接用的东西。1. 装好 Stata 之后先做这三件小事1.1 版本与安装包BE、SE、MP 到底差在哪很多人纠结的第一个问题是装哪个版本。这里要分清两个概念一个是版本号15、16、17、18 这样的年份号它决定语法特性和新命令另一个是 Edition也就是 BE、SE、MP它决定算力上限。语法上 BE 和 MP 基本一样差别在于能吃多大的数据。版本变量数上限矩阵维度上限并行计算适用场景Stata/BE798800无教学、小样本问卷Stata/SE3276711000无绝大多数科研项目Stata/MP6553411000多核千万行级数据、反复重抽样我自己的经验是如果你的数据只有几千行、几十个变量BE 完全够用跑回归的速度感受不到差别但一旦开始做上百万行的事务数据、或者要跑几百次 bootstrapMP 的多核优势非常明显。安装包从官方渠道或者所在机构的授权渠道获取。这一点我态度比较坚决不要用来路不明的压缩包。一是里面有捆绑脚本的风险二是版本混乱别人给的 do 文件在你这儿报一堆莫名其妙的错排查半天最后发现是版本不一致非常消耗精力。安装时注意把安装路径设成纯英文比如D:\Stata18不要放在我的文档\新建文件夹这种带中文和空格的路径下后面调用外部程序时会省掉很多麻烦。1.2 工作目录、日志与 do 文件让结果能复现我见过太多人是这样的在命令窗口敲一行、看一眼结果、再敲下一行第二天想复现完全记不起来昨天到底跑过什么。多元回归这种要反复试模型设定的活儿必须有个可复现的载体那就是 do 文件。新建一个项目文件夹结构大概长这样project/ data/ 原始数据只读不改 code/ do 文件 logs/ 运行日志 tables/ 导出的结果表 figures/ 图形输出每个 do 文件开头我会固定写这么几行version 17 clear all set more off set linesize 120 cd D:/project capture log close log using logs/20240115_main.log, replace text逐条说下为什么。version 17让脚本在更高版本下也按 17 的语法执行避免某个命令的行为在新版本里变了导致结果对不上。set more off是防止结果刷屏时卡在--more--等你按键批量跑的时候如果忘了这句脚本会永远停在那里。cd设置工作目录这样后面所有相对路径都不用写全。capture log close是防止上一次的日志没关干净导致log using报错capture会吞掉这个错误。log using ... , text里的text选项很重要它输出纯文本日志用任何编辑器都能打开比默认的 SMCL 格式通用得多。结尾别忘了log close1.3 中文变量名与编码的那些坑Stata 从 14 版开始全面支持 Unicode中文变量名在技术上是可以用的但我强烈建议变量名保持英文中文只放在变量标签label里。原因是变量名会出现在 do 文件、日志、导出表格的各个角落一旦跨系统Windows 到 Linux 服务器或者跨版本传递中文变量名出问题的概率远高于标签。数据导入常见的就是 Excel 和 CSVimport excel data/raw_2024.xlsx, sheet(Sheet1) firstrow clear import delimited data/raw_2024.csv, encoding(UTF-8) clearExcel 文件里那些男/女、是/否的文本列导入后是字符串型不能直接进回归需要转成数值encode gender, gen(gender_num) tab gender_numencode会按字符串的字母顺序自动编码并保留原来的文字作为值标签这样回归输出里显示的仍然是男/女而不是 1/2读表的时候舒服很多。如果是1、2这种被当成文本存进来的数字用destringdestring income, replace ignore(,)ignore(,)用来处理像 12,000 这种带千分位的写法。注意encode和destring一定要在合并数据之前做完而且要检查编码后的类别数对不对。我踩过的一次坑是问卷里缺失被录成了99直接进回归当成了一个巨大的合法数值把整个系数拉偏了。2. 多元回归到底在算什么先把直觉建立起来2.1 从控制变量这四个字说起假设你想比较不同城市的房价水平直接拿城市当自变量跑回归得出某城市房价显著更高这个结论基本没有意义——因为各地的人均收入、房屋面积、房龄结构都不一样你比的根本不是同一个东西。控制变量的作用就是把这些干扰因素拉到同一条水平线上再比。用数学说普通最小二乘OLS在做的是一件事找一组系数让预测值和真实值之间误差的平方和最小。min Σ(yᵢ − β₀ − β₁x₁ᵢ − … − βₖxₖᵢ)²为什么是平方而不是绝对值平方对大误差惩罚更重而且在数学上有解析解能直接写成矩阵公式一次算出来不需要迭代。代价是它对外侧的极端值特别敏感这也解释了后面为什么要专门处理异常值。系数怎么读以 β₁ 为例**在其他自变量保持不变的前提下x₁ 每增加一个单位y 平均变化 β₁ 个单位。**这句话里保持不变四个字是整个多元回归的灵魂报告结果时必须体现出来否则读者不知道你在控制什么。2.2 一张 regress 输出表逐块拆开看用 Stata 自带的数据跑一个最小例子sysuse auto, clear regress price mpg weight foreign输出分三块新手最容易只看中间那块其实上面和下面都藏着判断模型好坏的信息。第一块是方差分解表名称含义怎么看SS平方和Model Residual Totaldf自由度回归自由度 自变量个数MS均方SS / dfNumber of obs实际进入回归的样本量和原始数据行数不一致要查缺失F(k, n−k−1)整体显著性检验Prob F 小于 0.05 说明模型整体有效R-squared解释的方差比例横截面数据 0.2 以上就不算差Adj R-squared调整后 R²加变量后掉下去说明加的变量没用Root MSE残差标准差和 y 的量纲一致可理解为平均预测误差第二块是系数表每一行一个自变量五列分别是系数、标准误、t 值、p 值和 95% 置信区间。t 值就是系数除以标准误p 值是在真实系数为零的假设下观察到这么大 t 值的概率。第三块现在是空的但如果你在命令后加, vce(cluster id)之类的选项Stata 会把聚类信息、聚类数打在这里所以养成一个习惯跑完先往下翻一眼确认标准误类型是不是你以为的那种。2.3 显著性不是全部系数大小和置信区间有一个特别常见的误读p 小于 0.05 就万事大吉。实际上样本量够大时一个几乎可以忽略的微小效应也能显著。我见过一份十几万样本的数据某个变量系数只有 0.0003p 值小于 0.001作者兴奋地写显著正相关但换算成实际意义就是自变量翻一倍、因变量才动千分之几毫无业务价值。所以我现在的习惯是**先看系数大小和置信区间的宽度再看 p 值。**置信区间如果跨了很宽的范围说明估计得很不确定即使显著也不能下强结论。如果要把不同量纲的变量放在一起比重要性用标准化系数regress price mpg weight foreign, beta输出里多出一列 Beta它表示自变量每变化一个标准差、因变量变化多少个标准差这才是可比的影响力。不过强调一句标准化系数在论文里只能作为辅助主表还是用原始系数因为前者没法直接做经济或临床上的解释。3. 建模前的数据体检三张表不能省3.1 描述性统计与缺失值先搞清楚有多少样本能进模型很多人拿到数据直接regress跑完发现样本量从 5000 掉到 1800然后开始慌。多元回归是列表删除listwise deletion任何一个变量缺失这一行就被整体丢掉。所以建模前必须先看清楚缺失的分布。summarize price mpg weight foreign, detail tabstat price mpg weight foreign, stat(n mean sd min p50 max) col(stat) misstable summarize price mpg weight foreignsummarize加detail会给出百分位数、偏度、峰度快速判断分布形状。tabstat的好处是输出紧凑可以直接复制进论文的描述统计表。misstable summarize会显示每个变量的缺失数量如果某个变量缺失比例超过 10%就要认真考虑处理方案是删除这个变量、还是用多重插补mi系列命令。我个人的经验阈值是缺失低于 5% 且模式随机直接列表删除问题不大超过 10%最好做敏感性分析把插补结果和删除结果都跑一遍看看系数方向有没有变。提示描述统计和最终回归的样本必须一致。做法是在描述统计前先建立一个分析样本标记变量后面所有表都用它筛选避免出现描述统计用 5000 人、回归用 1800 人、审稿人一看就问你为什么不一致。3.2 相关矩阵与共线性初筛跑回归之前先看相关矩阵成本极低但收益很高pwcorr price mpg weight foreign, sig star(0.05) corr price mpg weight foreignpwcorr是两两相关corr是列表删除后的相关样本量不同的情况下两者结果会有差异一般汇报pwcorr并注明。看什么两件事。第一自变量和因变量的相关性方向是否符合预期如果某个变量的相关方向和你后面的回归系数方向相反这就是一个很好的警示信号通常意味着存在抑制效应或者共线性。第二自变量之间的相关系数有没有超过 0.8。超过这个量级两个变量基本在讲同一件事一起放进模型会导致标准误膨胀系数估计变得非常不稳定。我遇到过高相关变量时的一般处理顺序是先看业务含义哪个更贴近研究问题然后用 VIF 定量确认最后决定是删掉一个、还是合成一个综合指标、还是保留但在文中说明。绝不要因为删了可惜就把两个高度相关的变量硬塞进同一个模型。3.3 最大值、最小值与异常值几个必须会的命令这是实操里用得最多、但教程里往往一笔带过的一块。Stata 里取最大值最小值有好几种写法用途完全不同容易搞混。summarize price display r(min) display r(max) egen pmax max(price) egen pmin min(price) egen pmax_grp max(price), by(foreign) egen pmin_grp min(price), by(foreign) egen row_max rowmax(mpg weight length) egen row_min rowmin(mpg weight length)summarize之后用r(min)、r(max)取的是整体极值适合做逻辑校验比如年龄最大值是 999说明有录入错误。egen max()不加by()是整个样本的最大值加上by()就变成组内最大值这个在构造分组变量时特别常用。rowmax()、rowmin()是按行比较多个变量做量表题目的极端响应筛选时很好用。想看具体是哪几条记录最大最小用一个外部命令更省事ssc install extremes, replace extremes price mpg weight, n(10)它会直接列出最大和最小的若干条完整记录比sort price再list in 1/10舒服得多。发现异常值之后处理方式有三种选用哪种必须在文中交代清楚处理方式做法代价删除drop if price 20000样本量减少可能引入选择偏误缩尾winsorize把超过第 99 百分位的值替换为第 99 百分位的值改变尾部数值但保留样本截尾trim把超过阈值的观测直接删掉介于两者之间变量变换取对数、开方改变系数的解释方式缩尾最常用的实现ssc install winsor2, replace winsor2 price mpg, cuts(1 99) replace它同时生成新变量还是替换原变量由replace决定我建议不加replace生成price_w这样的新变量方便做敏感性对比。注意缩尾的比例不要凭感觉定。1% 和 99% 是常见默认值但如果你的数据本身尾部信息很重要比如收入分布、罕见病费用盲目缩尾会把真实效应抹掉。稳妥做法是分别用不缩尾、1/99 缩尾、5/95 缩尾跑三遍回归看系数方向是否一致。4. 把模型跑起来regress 及其实战变体4.1 基础语法、虚拟变量与交互项的写法regress的基本形式是regress 因变量 自变量1 自变量2 ...。真正让人头疼的是因子变量的写法记住四个符号就够用了regress price i.foreign c.mpg regress price c.mpg##c.weight regress price c.mpg##i.foreign regress price i.rep78##i.foreigni.表示把变量当分类变量处理自动生成虚拟变量c.表示当连续变量处理#表示只放交互项##表示既放主效应也放交互项。c.mpg##i.foreign这一行等价于同时放入了mpg、foreign和两者乘积三个项Stata 会替你算好。基准组的选择会影响所有虚拟变量系数的解释默认是取值最小的那一组。想改成别的regress price ib3.rep78##i.foreignib3.表示把 rep78 等于 3 的那组设为基准。这一步看着小但影响很大——基准组应该是样本量足够、含义最清晰的那一组否则所有系数都在跟一个奇怪的参照物比较。4.2 稳健标准误与聚类标准误什么时候必须加这是我在项目里最坚持的一件事默认永远先跑一遍robust再决定最终用哪种标准误。regress price mpg weight, robust regress price mpg weight, vce(cluster rep78)robust是异方差稳健标准误它不改变系数点估计只改变标准误和 p 值。也就是说加不加robust系数值是一个显著性可能完全不同。vce(cluster)是聚类稳健标准误用于观测值在某个维度上内部相关的场景比如同一家医院的多位患者、同一家公司的多个年份、同一个班级的多个学生。怎么判断该用哪个我的流程是这样先跑estat hettest做异方差检验。检验 p 值小于 0.05说明存在异方差用robust。数据有天然的嵌套结构个体嵌套于群体不管检验结果如何用vce(cluster 群体id)。提示聚类数太少时聚类标准误本身会严重低估真实变异。经验上聚类数少于 30 到 50 个时结果要谨慎对待可以考虑用 wild cluster bootstrap 做交叉验证。这条在审稿意见里经常被提提前做一遍能省一轮返修。4.3 面板数据与高维固定效应只要数据结构是同一个体多年或同一地区多期就该考虑固定效应它能把不随时间变化的个体异质性全部吸收掉这是横截面回归做不到的。xtset id year xtreg y x1 x2 i.year, fe vce(cluster id) areg y x1 x2 i.year, absorb(id) vce(cluster id) ssc install reghdfe, replace ssc install ftools, replace reghdfe y x1 x2, absorb(id year) vce(cluster id)xtreg, fe是最标准的写法。areg的好处是会报告组内固定效应的联合显著性而且对于只有一维固定效应的场景速度更快。reghdfe是我现在的主力因为它能一次吸收多个维度的固定效应比如同时吸收个体、年份、行业而且对内存的处理更聪明。选固定效应还是随机效应做 Hausman 检验。xtreg y x1 x2, fe estimates store fe_m xtreg y x1 x2, re estimates store re_m hausman fe_m re_m不过说实话在实际项目中我很少真的靠这个检验定结论。逻辑上如果你的研究问题是某个时变变量的效应而个体层面存在明显不可观测的异质性固定效应才是更安全的选择——它牺牲的是所有不随时间变化的变量的识别能力换来的是更强的因果解释力。代价要提前想清楚别跑完才发现核心变量因为不随时间变化被吸收掉了。4.4 结果导出别手敲系数手敲系数进表格是效率杀手而且极易出错。用esttabssc install estout, replace regress price mpg weight, robust estimates store m1 regress price mpg weight foreign, robust estimates store m2 regress price mpg weight foreign i.rep78, robust estimates store m3 esttab m1 m2 m3 using tables/reg_main.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a, labels(样本量 调整R2)) /// mtitles(模型1 模型2 模型3) nogapb(%9.3f)指定系数保留三位小数se()指定标准误star()自定义显著性星号的标准把 0.10 也标出来是很常见的做法。stats()决定底部报告哪些统计量。nogap去掉行间距表格更紧凑。如果要用 Word 的模板上报putexcel也可以但它需要你自己算好行列位置灵活性高、代价是代码量大。在项目周期紧张的时候我一般用esttab出 RTF再手工做最后的排版调整。5. 跑完不能直接交六项诊断与补救5.1 多重共线性VIF 只是起点regress price mpg weight length estat vif输出里每行一个变量VIF 和 1/VIF 两列。经验阈值是 VIF 大于 10对应 1/VIF 小于 0.1需要关注。但要注意VIF 只抓得到线性共线而且交互项天然就会带来很高的 VIF这并不是问题因为交互项和其构成变量本来就相关。我的处理方式对参与交互的连续变量先做中心化减去均值再构造交互项VIF 会好看很多而且系数解释更清晰——此时主效应表示的是在另一个变量取均值时的效应。summarize mpg generate mpg_c mpg - r(mean) regress price c.mpg_c##c.weight5.2 异方差与自相关别只做检验不管后果estat hettest estat imtest, whiteestat hettest是 Breusch-Pagan / Cook-Weisberg 检验estat imtest是 White 检验后者不假设异方差的具体形式更通用但自由度损失大小样本下功效偏低。发现异方差之后的处理优先级第一是检查函数形式是不是错了比如该取对数的没取很多时候异方差是模型设定问题的一个症状第二是用稳健标准误第三才是考虑加权最小二乘。面板数据还要查自相关ssc install xtserial, replace xtserial y x1 x2存在自相关时vce(cluster id)通常能同时缓解异方差和自相关问题这也是面板回归里聚类标准误几乎成了默认选项的原因。5.3 非线性、异常点与影响点回归诊断图比一堆检验统计量直观得多rvfplot lvr2plot predict resid, residuals predict cooksd, cooksd predict lev, leverage summarize cooksd, detailrvfplot是残差对拟合值图理想情况下应该是一团没有形状的随机散点。如果呈喇叭形提示异方差如果呈弯曲的 U 形说明遗漏了非线性项。lvr2plot是杠杆—残差平方图用来定位同时具有高杠杆和高残差的点。cooksd是 Cook 距离衡量单个观测对全部系数的影响经验阈值是 4/n超过的点值得逐条翻看原始记录。注意识别出影响点不等于可以删除。正确做法是把含该点和不含该点两个模型都跑一遍在文中报告两者说明核心结论是否稳健。直接删掉不报告是审稿时非常容易被抓的问题。5.4 内生性工具变量的最小可用流程如果核心自变量和误差项相关遗漏变量、反向因果、测量误差OLS 估计就是有偏的这时候要考虑工具变量。ivregress 2sls y x2 (x1 z1 z2), robust estat firststage estat endogenous estat overid三个检验看三件事estat firststage报告第一阶段的 F 值经验上小于 10 就要担心弱工具变量estat endogenous检验内生性是否真的存在如果不显著说明可以直接用 OLS没必要上工具变量estat overid是过度识别检验只有工具变量个数多于内生变量个数时才能做用来检验工具变量的外生性。我必须说清楚工具变量这套东西最大的难点从来不是命令而是有没有一个既满足相关性又满足外生性的工具。命令五分钟就能学会找到合适的工具变量可能要几个月。如果实在找不到退而求其次的做法是做敏感性分析比如用 Oster 检验评估遗漏变量需要多强才能推翻结论。6. 亚组分析分组回归还是交互项6.1 分组回归的三个陷阱最常见的亚组分析做法是按某个变量切分数据各跑一个回归regress y x1 x2 if gender 1, robust regress y x1 x2 if gender 0, robust然后看到男性组显著、女性组不显著就下结论说效应只存在于男性。这个推断在统计上是站不住的有三个问题第一两组样本量不同检验效能不同。一组显著、一组不显著很可能只是因为一组样本太少。第二你比较的是各自的 p 值而不是两组系数之间是否存在差异。第三做了多个亚组就产生了多重比较问题不校正的话假阳性率会明显上升。6.2 交互项加 margins我更推荐的做法正确做法是把亚组变量放进模型做交互regress y c.x1##i.gender x2 x3, robust contrast i.gender#c.x1 margins gender, dydx(x1) marginsplotcontrast i.gender#c.x1直接检验两组斜率是否相等这才是亚组差异的正式检验。margins gender, dydx(x1)给出每组的边际效应marginsplot把它画成图。要注意的是在regress里加交互项后margins的默认计算基准是其他变量取均值。如果你的模型里还有别的分类变量最好显式指定margins gender, dydx(x1) at(x2(10 20 30)) marginsplot, xdimension(x2)这样能得到边际效应如何随第三个变量变化的图用来展示效应异质性非常有说服力。6.3 做一张能看的亚组图如果亚组比较多比如三四个分组变量、每个变量两三水平一张汇总的森林图比一堆表格清晰得多。做法有两种一是用coefplot抽取系数和置信区间拼接二是把各亚组的margins结果导出到数据集手动构造点估计和上下限再用twoway rcap加scatter画出来。ssc install coefplot, replace regress y c.x1##i.gender##i.agegrp x2, robust coefplot, keep(*.gender#*.x1) xlabel(, labsize(small)) /// yline(0, lpattern(dash)) horizontalhorizontal让图横向排列yline(0)加一条零线是这类图的标配。多层交互会生成大量系数项建议不要一次画完按分组变量分批画。6.4 别把亚组分析和 Meta 分析混成一锅这里必须区分清楚两件经常被混淆的事。亚组分析处理的是个体层面数据用的是回归框架而 Meta 分析处理的是研究层面的效应量用的是加权平均框架。Stata 里的metan、network、mvmeta这些命令属于后者输入的是每个研究各自的效应量和标准误不是原始的病人数据。有人拿着只做了一次研究的个体数据硬套 Meta 分析的命令去跑这是方法论上的错位。正确的路径是如果你有多个独立研究的结果需要合并用 Meta 分析如果你有一份个体数据想探索效应在不同人群中的差异用交互项做亚组分析。两条路线各有各的假设和报告规范混着用很容易出问题。7. 报错与结果异常的排查清单7.1 常见报错速查Stata 的报错信息格式是r(错误码)新手看到一串数字就发懵。下面这张表覆盖了我日常碰到的绝大部分情况。错误码 / 提示常见原因处理方式r(111) variable not found变量名拼错、大小写不一致describe看实际变量名r(2000) no observations筛选条件把所有观测都排除了检查if条件和缺失值x omitted because of collinearity与已进入模型的变量完全共线检查是否重复放入、虚拟变量陷阱r(603) file not found路径写错、工作目录不对pwd看当前目录用绝对路径试试r(198) invalid syntax选项拼写错误、逗号位置不对看逗号前后选项之间用空格分隔type mismatch字符串变量进了数值运算destring或encoder(459) not sorted用by前没排序先sort或直接用bysortmissing values found变量有缺失但命令不允许加, missing选项或先处理缺失我自己的排查顺序是固定的**先describe确认变量在不在、什么类型再tab或summarize确认有没有数据最后才看命令语法。**这个顺序能覆盖八成以上的报错比盯着错误码猜快得多。7.2 结果不对劲的排查顺序报错好办最怕的是不报错但结果诡异。我的排查清单是这样的第一看样本量。跑完回归的Number of obs和你预期的差多少如果有大幅缩水一定是某个变量有缺失。第二看符号。核心变量的系数方向和你从文献、从常识得到的预期是否一致如果在描述性统计里是正相关在回归里变成负的大概率是共线性或者抑制效应。第三看量级。系数大得离谱比如因变量取值在 0 到 1 之间系数却是 300通常是单位问题——自变量有的是元、有的是万元混在一起了。第四看虚拟变量。某组的系数被 omitted检查是不是基准组设置不对或者某个水平在样本里完全没出现。第五换标准误再跑一遍。如果robust前后显著性天翻地覆说明存在异方差这本身就是一个需要报告的信息。7.3 数据量大到一定程度试试 ftools当数据从几万行涨到几百万行甚至千万行merge、collapse、duplicates这些操作的耗时会变得难以忍受。这时候用ftools包会快很多它是用 C 语言重写的一套替代命令内存占用和速度都有明显改善。ssc install ftools, replace fcollapse (mean) y (max) x1 (count) x2, by(id year) fmerge 1:1 id using data/big_aux.dta, nogen fegen tag tag(id) fegen n_obs count(y), by(id) fisid id year ftabulate gender agegrpfcollapse对应collapsefmerge对应mergefegen对应egenfisid用来快速检查某个键是不是唯一ftabulate对应tabulate。用法和原命令几乎一致改个前缀就能用学习成本很低。不过有一点要注意ftools的某些命令在极端情况下对缺失值和值标签的处理与官方命令略有差异。所以我会在正式跑之前先用小样本同时跑官方命令和ftools命令对比结果是否完全一致确认后再全量运行。这个交叉验证的习惯帮我抓到过好几次不对齐的问题尤其是涉及字符串键的合并。8. 一套我自己在用的回归项目模板8.1 目录结构与文件命名前面提过目录结构这里把命名规范也补上。我要求自己遵守两条do 文件按执行顺序编号输出文件带日期或者版本号。project/ 00_admin/ 说明文档、变量字典 01_data/ raw 只读clean 由代码生成 02_code/ 01_import.do 02_clean.do 03_desc.do 04_model.do 05_export.do 03_logs/ 自动生成不手改 04_output/ 回归表格、图形为什么坚持把 clean 数据也当成代码的产物而不是手工修改的结果因为一旦你手动改过某个数据文件后面再想复现就彻底断了。所有清洗逻辑都必须写在代码里哪怕只是一次性修改一个错值也写成replace x 30 if id 1024这样别人包括半年后的你自己才能看懂发生了什么。8.2 do 文件骨架主模型文件大概长这样我基本每个项目都是这个套路version 17 clear all set more off set seed 20240115 cd D:/project capture log close log using 03_logs/04_model.log, replace text use 01_data/clean_2024.dta, clear global xvars x1 x2 x3 x4 global controls age gender edu income regress y $xvars $controls, robust estimates store m1 regress y $xvars $controls i.year i.region, robust estimates store m2 regress y c.x1##i.gender $controls, robust estimates store m3 esttab m1 m2 m3 using 04_output/reg_main.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a, labels(样本量 调整R2)) /// mtitles(基准 加固定效应 交互) log close用global宏来管理变量列表是我做久了之后最大的效率提升点。变量清单在文件开头集中定义后面所有模型都引用$controls需要改控制变量时只改一行不用在几十个模型里逐个替换。set seed放在开头看着多余但只要后面有任何涉及随机数的操作bootstrap、多重插补、随机抽样没有固定种子就意味着结果不可复现。养成习惯比事后补救容易得多。8.3 参数化与批量跑模型如果需要跑几十个模型比如对十几个因变量分别做同样的回归用循环比复制粘贴靠谱foreach y in y1 y2 y3 y4 { regress y $xvars $controls, robust estimates store m_y } esttab m_* using 04_output/reg_batch.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a) mtitlesforeach里的y是局部宏引用注意是反引号加单引号写错一个符号就会报invalid syntax这是新手最常见的错误之一。还有一个我特别推荐的技巧跑完所有模型后把关键系数统一导出来做趋势检查。tempname fh file open fh using 04_output/coef_trend.csv, write replace file write fh model,var,coef,se,p _n foreach m in m1 m2 m3 { estimates restore m foreach v in x1 x2 { file write fh m,v, /// %9.4f (_b[v]) , %9.4f (_se[v]) , /// %9.4f (2*ttail(e(df_r), abs(_b[v]/_se[v]))) _n } } file close fh这样得到的是一张模型—变量—系数—标准误—p 值的长表用任何画图工具都能看出核心变量系数在不同设定下是否稳定。系数稳定性检查比单个模型的各种检验都更能说明结果靠不靠谱这个习惯我是从一个做应用微观的同行那儿学来的之后所有项目都保留了下来。关于_b[]和_se[]它们只能在估计完成后立即访问所以必须配合estimates restore使用。另外注意e(df_r)是残差自由度只有在用robust或cluster选项时它才存在普通 OLS 下要换成e(N) - e(df_m) - 1自己算。这个小细节不注意的话p 值会全变成缺失。最后再分享一个我在实操中体会比较深的点**多元回归的功夫八成花在跑回归之前。**模型的设定、变量的取舍、样本的选择、缺失的处理这些决定远比用什么命令、加不加某个选项重要。我早期也走过一段追命令的弯路到处找最新最炫的估计方法后来才明白一个设定干净、诊断完整、结果稳健的普通 OLS比一个堆满高级命令但连描述统计都对不上的模型可信度高得多。真要把这套流程跑顺建议拿一份自己熟悉的数据从描述统计开始完整走一遍把日志留下来过两周再回头看那段代码能看懂、能复现基本功就算过关了。
Stata多元回归分析实操：数据清洗、稳健标准误与亚组交互
很多人第一次接触 Stata是因为导师或同事丢过来一个 do 文件说照着跑一遍就行。我自己当年也是这么入门的一份几十个变量的问卷数据一张表格要填十几个模型的系数用表格软件硬算根本算不出来。后来项目做得多了从横截面问卷到面板数据从只盯 p 值到认真处理内生性和异质性多元回归分析在Stata里几乎成了我每周都要碰的东西。这篇东西不打算复述教科书只讲实操里踩过、绕过、最后固定下来的一整套流程装完之后先做什么、regress输出那张表每一列到底在说什么、稳健标准误什么时候必须加、亚组分析该分组回归还是做交互、报错信息怎么读以及一堆看起来不起眼但能省半天时间的细节命令。刚上手的学生能照着走完一遍做了几年想回头补基础的从业者也能从里面挑到能直接用的东西。1. 装好 Stata 之后先做这三件小事1.1 版本与安装包BE、SE、MP 到底差在哪很多人纠结的第一个问题是装哪个版本。这里要分清两个概念一个是版本号15、16、17、18 这样的年份号它决定语法特性和新命令另一个是 Edition也就是 BE、SE、MP它决定算力上限。语法上 BE 和 MP 基本一样差别在于能吃多大的数据。版本变量数上限矩阵维度上限并行计算适用场景Stata/BE798800无教学、小样本问卷Stata/SE3276711000无绝大多数科研项目Stata/MP6553411000多核千万行级数据、反复重抽样我自己的经验是如果你的数据只有几千行、几十个变量BE 完全够用跑回归的速度感受不到差别但一旦开始做上百万行的事务数据、或者要跑几百次 bootstrapMP 的多核优势非常明显。安装包从官方渠道或者所在机构的授权渠道获取。这一点我态度比较坚决不要用来路不明的压缩包。一是里面有捆绑脚本的风险二是版本混乱别人给的 do 文件在你这儿报一堆莫名其妙的错排查半天最后发现是版本不一致非常消耗精力。安装时注意把安装路径设成纯英文比如D:\Stata18不要放在我的文档\新建文件夹这种带中文和空格的路径下后面调用外部程序时会省掉很多麻烦。1.2 工作目录、日志与 do 文件让结果能复现我见过太多人是这样的在命令窗口敲一行、看一眼结果、再敲下一行第二天想复现完全记不起来昨天到底跑过什么。多元回归这种要反复试模型设定的活儿必须有个可复现的载体那就是 do 文件。新建一个项目文件夹结构大概长这样project/ data/ 原始数据只读不改 code/ do 文件 logs/ 运行日志 tables/ 导出的结果表 figures/ 图形输出每个 do 文件开头我会固定写这么几行version 17 clear all set more off set linesize 120 cd D:/project capture log close log using logs/20240115_main.log, replace text逐条说下为什么。version 17让脚本在更高版本下也按 17 的语法执行避免某个命令的行为在新版本里变了导致结果对不上。set more off是防止结果刷屏时卡在--more--等你按键批量跑的时候如果忘了这句脚本会永远停在那里。cd设置工作目录这样后面所有相对路径都不用写全。capture log close是防止上一次的日志没关干净导致log using报错capture会吞掉这个错误。log using ... , text里的text选项很重要它输出纯文本日志用任何编辑器都能打开比默认的 SMCL 格式通用得多。结尾别忘了log close1.3 中文变量名与编码的那些坑Stata 从 14 版开始全面支持 Unicode中文变量名在技术上是可以用的但我强烈建议变量名保持英文中文只放在变量标签label里。原因是变量名会出现在 do 文件、日志、导出表格的各个角落一旦跨系统Windows 到 Linux 服务器或者跨版本传递中文变量名出问题的概率远高于标签。数据导入常见的就是 Excel 和 CSVimport excel data/raw_2024.xlsx, sheet(Sheet1) firstrow clear import delimited data/raw_2024.csv, encoding(UTF-8) clearExcel 文件里那些男/女、是/否的文本列导入后是字符串型不能直接进回归需要转成数值encode gender, gen(gender_num) tab gender_numencode会按字符串的字母顺序自动编码并保留原来的文字作为值标签这样回归输出里显示的仍然是男/女而不是 1/2读表的时候舒服很多。如果是1、2这种被当成文本存进来的数字用destringdestring income, replace ignore(,)ignore(,)用来处理像 12,000 这种带千分位的写法。注意encode和destring一定要在合并数据之前做完而且要检查编码后的类别数对不对。我踩过的一次坑是问卷里缺失被录成了99直接进回归当成了一个巨大的合法数值把整个系数拉偏了。2. 多元回归到底在算什么先把直觉建立起来2.1 从控制变量这四个字说起假设你想比较不同城市的房价水平直接拿城市当自变量跑回归得出某城市房价显著更高这个结论基本没有意义——因为各地的人均收入、房屋面积、房龄结构都不一样你比的根本不是同一个东西。控制变量的作用就是把这些干扰因素拉到同一条水平线上再比。用数学说普通最小二乘OLS在做的是一件事找一组系数让预测值和真实值之间误差的平方和最小。min Σ(yᵢ − β₀ − β₁x₁ᵢ − … − βₖxₖᵢ)²为什么是平方而不是绝对值平方对大误差惩罚更重而且在数学上有解析解能直接写成矩阵公式一次算出来不需要迭代。代价是它对外侧的极端值特别敏感这也解释了后面为什么要专门处理异常值。系数怎么读以 β₁ 为例**在其他自变量保持不变的前提下x₁ 每增加一个单位y 平均变化 β₁ 个单位。**这句话里保持不变四个字是整个多元回归的灵魂报告结果时必须体现出来否则读者不知道你在控制什么。2.2 一张 regress 输出表逐块拆开看用 Stata 自带的数据跑一个最小例子sysuse auto, clear regress price mpg weight foreign输出分三块新手最容易只看中间那块其实上面和下面都藏着判断模型好坏的信息。第一块是方差分解表名称含义怎么看SS平方和Model Residual Totaldf自由度回归自由度 自变量个数MS均方SS / dfNumber of obs实际进入回归的样本量和原始数据行数不一致要查缺失F(k, n−k−1)整体显著性检验Prob F 小于 0.05 说明模型整体有效R-squared解释的方差比例横截面数据 0.2 以上就不算差Adj R-squared调整后 R²加变量后掉下去说明加的变量没用Root MSE残差标准差和 y 的量纲一致可理解为平均预测误差第二块是系数表每一行一个自变量五列分别是系数、标准误、t 值、p 值和 95% 置信区间。t 值就是系数除以标准误p 值是在真实系数为零的假设下观察到这么大 t 值的概率。第三块现在是空的但如果你在命令后加, vce(cluster id)之类的选项Stata 会把聚类信息、聚类数打在这里所以养成一个习惯跑完先往下翻一眼确认标准误类型是不是你以为的那种。2.3 显著性不是全部系数大小和置信区间有一个特别常见的误读p 小于 0.05 就万事大吉。实际上样本量够大时一个几乎可以忽略的微小效应也能显著。我见过一份十几万样本的数据某个变量系数只有 0.0003p 值小于 0.001作者兴奋地写显著正相关但换算成实际意义就是自变量翻一倍、因变量才动千分之几毫无业务价值。所以我现在的习惯是**先看系数大小和置信区间的宽度再看 p 值。**置信区间如果跨了很宽的范围说明估计得很不确定即使显著也不能下强结论。如果要把不同量纲的变量放在一起比重要性用标准化系数regress price mpg weight foreign, beta输出里多出一列 Beta它表示自变量每变化一个标准差、因变量变化多少个标准差这才是可比的影响力。不过强调一句标准化系数在论文里只能作为辅助主表还是用原始系数因为前者没法直接做经济或临床上的解释。3. 建模前的数据体检三张表不能省3.1 描述性统计与缺失值先搞清楚有多少样本能进模型很多人拿到数据直接regress跑完发现样本量从 5000 掉到 1800然后开始慌。多元回归是列表删除listwise deletion任何一个变量缺失这一行就被整体丢掉。所以建模前必须先看清楚缺失的分布。summarize price mpg weight foreign, detail tabstat price mpg weight foreign, stat(n mean sd min p50 max) col(stat) misstable summarize price mpg weight foreignsummarize加detail会给出百分位数、偏度、峰度快速判断分布形状。tabstat的好处是输出紧凑可以直接复制进论文的描述统计表。misstable summarize会显示每个变量的缺失数量如果某个变量缺失比例超过 10%就要认真考虑处理方案是删除这个变量、还是用多重插补mi系列命令。我个人的经验阈值是缺失低于 5% 且模式随机直接列表删除问题不大超过 10%最好做敏感性分析把插补结果和删除结果都跑一遍看看系数方向有没有变。提示描述统计和最终回归的样本必须一致。做法是在描述统计前先建立一个分析样本标记变量后面所有表都用它筛选避免出现描述统计用 5000 人、回归用 1800 人、审稿人一看就问你为什么不一致。3.2 相关矩阵与共线性初筛跑回归之前先看相关矩阵成本极低但收益很高pwcorr price mpg weight foreign, sig star(0.05) corr price mpg weight foreignpwcorr是两两相关corr是列表删除后的相关样本量不同的情况下两者结果会有差异一般汇报pwcorr并注明。看什么两件事。第一自变量和因变量的相关性方向是否符合预期如果某个变量的相关方向和你后面的回归系数方向相反这就是一个很好的警示信号通常意味着存在抑制效应或者共线性。第二自变量之间的相关系数有没有超过 0.8。超过这个量级两个变量基本在讲同一件事一起放进模型会导致标准误膨胀系数估计变得非常不稳定。我遇到过高相关变量时的一般处理顺序是先看业务含义哪个更贴近研究问题然后用 VIF 定量确认最后决定是删掉一个、还是合成一个综合指标、还是保留但在文中说明。绝不要因为删了可惜就把两个高度相关的变量硬塞进同一个模型。3.3 最大值、最小值与异常值几个必须会的命令这是实操里用得最多、但教程里往往一笔带过的一块。Stata 里取最大值最小值有好几种写法用途完全不同容易搞混。summarize price display r(min) display r(max) egen pmax max(price) egen pmin min(price) egen pmax_grp max(price), by(foreign) egen pmin_grp min(price), by(foreign) egen row_max rowmax(mpg weight length) egen row_min rowmin(mpg weight length)summarize之后用r(min)、r(max)取的是整体极值适合做逻辑校验比如年龄最大值是 999说明有录入错误。egen max()不加by()是整个样本的最大值加上by()就变成组内最大值这个在构造分组变量时特别常用。rowmax()、rowmin()是按行比较多个变量做量表题目的极端响应筛选时很好用。想看具体是哪几条记录最大最小用一个外部命令更省事ssc install extremes, replace extremes price mpg weight, n(10)它会直接列出最大和最小的若干条完整记录比sort price再list in 1/10舒服得多。发现异常值之后处理方式有三种选用哪种必须在文中交代清楚处理方式做法代价删除drop if price 20000样本量减少可能引入选择偏误缩尾winsorize把超过第 99 百分位的值替换为第 99 百分位的值改变尾部数值但保留样本截尾trim把超过阈值的观测直接删掉介于两者之间变量变换取对数、开方改变系数的解释方式缩尾最常用的实现ssc install winsor2, replace winsor2 price mpg, cuts(1 99) replace它同时生成新变量还是替换原变量由replace决定我建议不加replace生成price_w这样的新变量方便做敏感性对比。注意缩尾的比例不要凭感觉定。1% 和 99% 是常见默认值但如果你的数据本身尾部信息很重要比如收入分布、罕见病费用盲目缩尾会把真实效应抹掉。稳妥做法是分别用不缩尾、1/99 缩尾、5/95 缩尾跑三遍回归看系数方向是否一致。4. 把模型跑起来regress 及其实战变体4.1 基础语法、虚拟变量与交互项的写法regress的基本形式是regress 因变量 自变量1 自变量2 ...。真正让人头疼的是因子变量的写法记住四个符号就够用了regress price i.foreign c.mpg regress price c.mpg##c.weight regress price c.mpg##i.foreign regress price i.rep78##i.foreigni.表示把变量当分类变量处理自动生成虚拟变量c.表示当连续变量处理#表示只放交互项##表示既放主效应也放交互项。c.mpg##i.foreign这一行等价于同时放入了mpg、foreign和两者乘积三个项Stata 会替你算好。基准组的选择会影响所有虚拟变量系数的解释默认是取值最小的那一组。想改成别的regress price ib3.rep78##i.foreignib3.表示把 rep78 等于 3 的那组设为基准。这一步看着小但影响很大——基准组应该是样本量足够、含义最清晰的那一组否则所有系数都在跟一个奇怪的参照物比较。4.2 稳健标准误与聚类标准误什么时候必须加这是我在项目里最坚持的一件事默认永远先跑一遍robust再决定最终用哪种标准误。regress price mpg weight, robust regress price mpg weight, vce(cluster rep78)robust是异方差稳健标准误它不改变系数点估计只改变标准误和 p 值。也就是说加不加robust系数值是一个显著性可能完全不同。vce(cluster)是聚类稳健标准误用于观测值在某个维度上内部相关的场景比如同一家医院的多位患者、同一家公司的多个年份、同一个班级的多个学生。怎么判断该用哪个我的流程是这样先跑estat hettest做异方差检验。检验 p 值小于 0.05说明存在异方差用robust。数据有天然的嵌套结构个体嵌套于群体不管检验结果如何用vce(cluster 群体id)。提示聚类数太少时聚类标准误本身会严重低估真实变异。经验上聚类数少于 30 到 50 个时结果要谨慎对待可以考虑用 wild cluster bootstrap 做交叉验证。这条在审稿意见里经常被提提前做一遍能省一轮返修。4.3 面板数据与高维固定效应只要数据结构是同一个体多年或同一地区多期就该考虑固定效应它能把不随时间变化的个体异质性全部吸收掉这是横截面回归做不到的。xtset id year xtreg y x1 x2 i.year, fe vce(cluster id) areg y x1 x2 i.year, absorb(id) vce(cluster id) ssc install reghdfe, replace ssc install ftools, replace reghdfe y x1 x2, absorb(id year) vce(cluster id)xtreg, fe是最标准的写法。areg的好处是会报告组内固定效应的联合显著性而且对于只有一维固定效应的场景速度更快。reghdfe是我现在的主力因为它能一次吸收多个维度的固定效应比如同时吸收个体、年份、行业而且对内存的处理更聪明。选固定效应还是随机效应做 Hausman 检验。xtreg y x1 x2, fe estimates store fe_m xtreg y x1 x2, re estimates store re_m hausman fe_m re_m不过说实话在实际项目中我很少真的靠这个检验定结论。逻辑上如果你的研究问题是某个时变变量的效应而个体层面存在明显不可观测的异质性固定效应才是更安全的选择——它牺牲的是所有不随时间变化的变量的识别能力换来的是更强的因果解释力。代价要提前想清楚别跑完才发现核心变量因为不随时间变化被吸收掉了。4.4 结果导出别手敲系数手敲系数进表格是效率杀手而且极易出错。用esttabssc install estout, replace regress price mpg weight, robust estimates store m1 regress price mpg weight foreign, robust estimates store m2 regress price mpg weight foreign i.rep78, robust estimates store m3 esttab m1 m2 m3 using tables/reg_main.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a, labels(样本量 调整R2)) /// mtitles(模型1 模型2 模型3) nogapb(%9.3f)指定系数保留三位小数se()指定标准误star()自定义显著性星号的标准把 0.10 也标出来是很常见的做法。stats()决定底部报告哪些统计量。nogap去掉行间距表格更紧凑。如果要用 Word 的模板上报putexcel也可以但它需要你自己算好行列位置灵活性高、代价是代码量大。在项目周期紧张的时候我一般用esttab出 RTF再手工做最后的排版调整。5. 跑完不能直接交六项诊断与补救5.1 多重共线性VIF 只是起点regress price mpg weight length estat vif输出里每行一个变量VIF 和 1/VIF 两列。经验阈值是 VIF 大于 10对应 1/VIF 小于 0.1需要关注。但要注意VIF 只抓得到线性共线而且交互项天然就会带来很高的 VIF这并不是问题因为交互项和其构成变量本来就相关。我的处理方式对参与交互的连续变量先做中心化减去均值再构造交互项VIF 会好看很多而且系数解释更清晰——此时主效应表示的是在另一个变量取均值时的效应。summarize mpg generate mpg_c mpg - r(mean) regress price c.mpg_c##c.weight5.2 异方差与自相关别只做检验不管后果estat hettest estat imtest, whiteestat hettest是 Breusch-Pagan / Cook-Weisberg 检验estat imtest是 White 检验后者不假设异方差的具体形式更通用但自由度损失大小样本下功效偏低。发现异方差之后的处理优先级第一是检查函数形式是不是错了比如该取对数的没取很多时候异方差是模型设定问题的一个症状第二是用稳健标准误第三才是考虑加权最小二乘。面板数据还要查自相关ssc install xtserial, replace xtserial y x1 x2存在自相关时vce(cluster id)通常能同时缓解异方差和自相关问题这也是面板回归里聚类标准误几乎成了默认选项的原因。5.3 非线性、异常点与影响点回归诊断图比一堆检验统计量直观得多rvfplot lvr2plot predict resid, residuals predict cooksd, cooksd predict lev, leverage summarize cooksd, detailrvfplot是残差对拟合值图理想情况下应该是一团没有形状的随机散点。如果呈喇叭形提示异方差如果呈弯曲的 U 形说明遗漏了非线性项。lvr2plot是杠杆—残差平方图用来定位同时具有高杠杆和高残差的点。cooksd是 Cook 距离衡量单个观测对全部系数的影响经验阈值是 4/n超过的点值得逐条翻看原始记录。注意识别出影响点不等于可以删除。正确做法是把含该点和不含该点两个模型都跑一遍在文中报告两者说明核心结论是否稳健。直接删掉不报告是审稿时非常容易被抓的问题。5.4 内生性工具变量的最小可用流程如果核心自变量和误差项相关遗漏变量、反向因果、测量误差OLS 估计就是有偏的这时候要考虑工具变量。ivregress 2sls y x2 (x1 z1 z2), robust estat firststage estat endogenous estat overid三个检验看三件事estat firststage报告第一阶段的 F 值经验上小于 10 就要担心弱工具变量estat endogenous检验内生性是否真的存在如果不显著说明可以直接用 OLS没必要上工具变量estat overid是过度识别检验只有工具变量个数多于内生变量个数时才能做用来检验工具变量的外生性。我必须说清楚工具变量这套东西最大的难点从来不是命令而是有没有一个既满足相关性又满足外生性的工具。命令五分钟就能学会找到合适的工具变量可能要几个月。如果实在找不到退而求其次的做法是做敏感性分析比如用 Oster 检验评估遗漏变量需要多强才能推翻结论。6. 亚组分析分组回归还是交互项6.1 分组回归的三个陷阱最常见的亚组分析做法是按某个变量切分数据各跑一个回归regress y x1 x2 if gender 1, robust regress y x1 x2 if gender 0, robust然后看到男性组显著、女性组不显著就下结论说效应只存在于男性。这个推断在统计上是站不住的有三个问题第一两组样本量不同检验效能不同。一组显著、一组不显著很可能只是因为一组样本太少。第二你比较的是各自的 p 值而不是两组系数之间是否存在差异。第三做了多个亚组就产生了多重比较问题不校正的话假阳性率会明显上升。6.2 交互项加 margins我更推荐的做法正确做法是把亚组变量放进模型做交互regress y c.x1##i.gender x2 x3, robust contrast i.gender#c.x1 margins gender, dydx(x1) marginsplotcontrast i.gender#c.x1直接检验两组斜率是否相等这才是亚组差异的正式检验。margins gender, dydx(x1)给出每组的边际效应marginsplot把它画成图。要注意的是在regress里加交互项后margins的默认计算基准是其他变量取均值。如果你的模型里还有别的分类变量最好显式指定margins gender, dydx(x1) at(x2(10 20 30)) marginsplot, xdimension(x2)这样能得到边际效应如何随第三个变量变化的图用来展示效应异质性非常有说服力。6.3 做一张能看的亚组图如果亚组比较多比如三四个分组变量、每个变量两三水平一张汇总的森林图比一堆表格清晰得多。做法有两种一是用coefplot抽取系数和置信区间拼接二是把各亚组的margins结果导出到数据集手动构造点估计和上下限再用twoway rcap加scatter画出来。ssc install coefplot, replace regress y c.x1##i.gender##i.agegrp x2, robust coefplot, keep(*.gender#*.x1) xlabel(, labsize(small)) /// yline(0, lpattern(dash)) horizontalhorizontal让图横向排列yline(0)加一条零线是这类图的标配。多层交互会生成大量系数项建议不要一次画完按分组变量分批画。6.4 别把亚组分析和 Meta 分析混成一锅这里必须区分清楚两件经常被混淆的事。亚组分析处理的是个体层面数据用的是回归框架而 Meta 分析处理的是研究层面的效应量用的是加权平均框架。Stata 里的metan、network、mvmeta这些命令属于后者输入的是每个研究各自的效应量和标准误不是原始的病人数据。有人拿着只做了一次研究的个体数据硬套 Meta 分析的命令去跑这是方法论上的错位。正确的路径是如果你有多个独立研究的结果需要合并用 Meta 分析如果你有一份个体数据想探索效应在不同人群中的差异用交互项做亚组分析。两条路线各有各的假设和报告规范混着用很容易出问题。7. 报错与结果异常的排查清单7.1 常见报错速查Stata 的报错信息格式是r(错误码)新手看到一串数字就发懵。下面这张表覆盖了我日常碰到的绝大部分情况。错误码 / 提示常见原因处理方式r(111) variable not found变量名拼错、大小写不一致describe看实际变量名r(2000) no observations筛选条件把所有观测都排除了检查if条件和缺失值x omitted because of collinearity与已进入模型的变量完全共线检查是否重复放入、虚拟变量陷阱r(603) file not found路径写错、工作目录不对pwd看当前目录用绝对路径试试r(198) invalid syntax选项拼写错误、逗号位置不对看逗号前后选项之间用空格分隔type mismatch字符串变量进了数值运算destring或encoder(459) not sorted用by前没排序先sort或直接用bysortmissing values found变量有缺失但命令不允许加, missing选项或先处理缺失我自己的排查顺序是固定的**先describe确认变量在不在、什么类型再tab或summarize确认有没有数据最后才看命令语法。**这个顺序能覆盖八成以上的报错比盯着错误码猜快得多。7.2 结果不对劲的排查顺序报错好办最怕的是不报错但结果诡异。我的排查清单是这样的第一看样本量。跑完回归的Number of obs和你预期的差多少如果有大幅缩水一定是某个变量有缺失。第二看符号。核心变量的系数方向和你从文献、从常识得到的预期是否一致如果在描述性统计里是正相关在回归里变成负的大概率是共线性或者抑制效应。第三看量级。系数大得离谱比如因变量取值在 0 到 1 之间系数却是 300通常是单位问题——自变量有的是元、有的是万元混在一起了。第四看虚拟变量。某组的系数被 omitted检查是不是基准组设置不对或者某个水平在样本里完全没出现。第五换标准误再跑一遍。如果robust前后显著性天翻地覆说明存在异方差这本身就是一个需要报告的信息。7.3 数据量大到一定程度试试 ftools当数据从几万行涨到几百万行甚至千万行merge、collapse、duplicates这些操作的耗时会变得难以忍受。这时候用ftools包会快很多它是用 C 语言重写的一套替代命令内存占用和速度都有明显改善。ssc install ftools, replace fcollapse (mean) y (max) x1 (count) x2, by(id year) fmerge 1:1 id using data/big_aux.dta, nogen fegen tag tag(id) fegen n_obs count(y), by(id) fisid id year ftabulate gender agegrpfcollapse对应collapsefmerge对应mergefegen对应egenfisid用来快速检查某个键是不是唯一ftabulate对应tabulate。用法和原命令几乎一致改个前缀就能用学习成本很低。不过有一点要注意ftools的某些命令在极端情况下对缺失值和值标签的处理与官方命令略有差异。所以我会在正式跑之前先用小样本同时跑官方命令和ftools命令对比结果是否完全一致确认后再全量运行。这个交叉验证的习惯帮我抓到过好几次不对齐的问题尤其是涉及字符串键的合并。8. 一套我自己在用的回归项目模板8.1 目录结构与文件命名前面提过目录结构这里把命名规范也补上。我要求自己遵守两条do 文件按执行顺序编号输出文件带日期或者版本号。project/ 00_admin/ 说明文档、变量字典 01_data/ raw 只读clean 由代码生成 02_code/ 01_import.do 02_clean.do 03_desc.do 04_model.do 05_export.do 03_logs/ 自动生成不手改 04_output/ 回归表格、图形为什么坚持把 clean 数据也当成代码的产物而不是手工修改的结果因为一旦你手动改过某个数据文件后面再想复现就彻底断了。所有清洗逻辑都必须写在代码里哪怕只是一次性修改一个错值也写成replace x 30 if id 1024这样别人包括半年后的你自己才能看懂发生了什么。8.2 do 文件骨架主模型文件大概长这样我基本每个项目都是这个套路version 17 clear all set more off set seed 20240115 cd D:/project capture log close log using 03_logs/04_model.log, replace text use 01_data/clean_2024.dta, clear global xvars x1 x2 x3 x4 global controls age gender edu income regress y $xvars $controls, robust estimates store m1 regress y $xvars $controls i.year i.region, robust estimates store m2 regress y c.x1##i.gender $controls, robust estimates store m3 esttab m1 m2 m3 using 04_output/reg_main.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a, labels(样本量 调整R2)) /// mtitles(基准 加固定效应 交互) log close用global宏来管理变量列表是我做久了之后最大的效率提升点。变量清单在文件开头集中定义后面所有模型都引用$controls需要改控制变量时只改一行不用在几十个模型里逐个替换。set seed放在开头看着多余但只要后面有任何涉及随机数的操作bootstrap、多重插补、随机抽样没有固定种子就意味着结果不可复现。养成习惯比事后补救容易得多。8.3 参数化与批量跑模型如果需要跑几十个模型比如对十几个因变量分别做同样的回归用循环比复制粘贴靠谱foreach y in y1 y2 y3 y4 { regress y $xvars $controls, robust estimates store m_y } esttab m_* using 04_output/reg_batch.rtf, replace /// b(%9.3f) se(%9.3f) star(* 0.10 ** 0.05 *** 0.01) /// stats(N r2_a) mtitlesforeach里的y是局部宏引用注意是反引号加单引号写错一个符号就会报invalid syntax这是新手最常见的错误之一。还有一个我特别推荐的技巧跑完所有模型后把关键系数统一导出来做趋势检查。tempname fh file open fh using 04_output/coef_trend.csv, write replace file write fh model,var,coef,se,p _n foreach m in m1 m2 m3 { estimates restore m foreach v in x1 x2 { file write fh m,v, /// %9.4f (_b[v]) , %9.4f (_se[v]) , /// %9.4f (2*ttail(e(df_r), abs(_b[v]/_se[v]))) _n } } file close fh这样得到的是一张模型—变量—系数—标准误—p 值的长表用任何画图工具都能看出核心变量系数在不同设定下是否稳定。系数稳定性检查比单个模型的各种检验都更能说明结果靠不靠谱这个习惯我是从一个做应用微观的同行那儿学来的之后所有项目都保留了下来。关于_b[]和_se[]它们只能在估计完成后立即访问所以必须配合estimates restore使用。另外注意e(df_r)是残差自由度只有在用robust或cluster选项时它才存在普通 OLS 下要换成e(N) - e(df_m) - 1自己算。这个小细节不注意的话p 值会全变成缺失。最后再分享一个我在实操中体会比较深的点**多元回归的功夫八成花在跑回归之前。**模型的设定、变量的取舍、样本的选择、缺失的处理这些决定远比用什么命令、加不加某个选项重要。我早期也走过一段追命令的弯路到处找最新最炫的估计方法后来才明白一个设定干净、诊断完整、结果稳健的普通 OLS比一个堆满高级命令但连描述统计都对不上的模型可信度高得多。真要把这套流程跑顺建议拿一份自己熟悉的数据从描述统计开始完整走一遍把日志留下来过两周再回头看那段代码能看懂、能复现基本功就算过关了。
