上次给大家带来了使用 esttab 命令从 Stata 导出 描述性统计表、相关系数矩阵、分组 t 检验表、多元多模型回归结果。并且详细教大家怎么定制化自己的表格，使其更符合学术论文三线表要求。详情可见
对于一般的朋友，这已经足以应对他们的要求了，但是，总有一些朋(chu)友(nv)们(zuo)会有更(bian)高(tai)的要求。本期将为大家带来更加详细的介绍。本期为大家带来多方程模型结果导出（如似不相关回归，分组系数检验），变量分组展示（如控制了个人、企业、行业层面的变量，希望分组展示），以及矩阵导出（bootstrap法检验中介效应结果导出）
先补一个坑
上一篇中，细心的小伙伴应该也发现了，尽管我们使用中文的变量 label 可以让导出的变量以中文形式出现，但是常数项 Constant 依然保持英文，我们该怎么让它也输出 常数项 这个中文呢？
其实，esttab 本身也已经提供了相关的方法，也就是 coeflabels() 选项，可以简写为 coef()，代码如下（coef()被我放在esttab的 note选项后面）
* --------------------------------------------------
* 与上一次相同的设定
* --------------------------------------------------
use "http://www.stata-press.com/data/r16/nlswork.dta",clear
xtset idcode year // 指定面板数据
glo y ln_wage // 因变量
glo x hours // 自变量
glo z c_city // 工具变量
glo m race // 调节变量
glo c grade age // 控制变量
glo f i.year i.ind_code // 年份和行业固定效应
// 变量标签标记为中文
label var $y "工资"
label var $x "工时"
label var $z "中心城市"
label var $m "种族"
label var grade "学历"
label var age "年龄"
// 值标签标记为中文
label define racelbl 1 "白人" 2 "黑人" 3 "其他", replace
* --------------------------------------------------
* 储存回归结果
* --------------------------------------------------
eststo olsbase: reg $y $c $f , vce(cluster idcode)
eststo olsmain: reg $y $x $c $f , vce(cluster idcode)
eststo olsint : reg $y c.${x}##ib3.${m} $c $f , vce(cluster idcode)
eststo regivm : ivregress 2sls $y (${x} = ${z}) $c $f , vce(cluster idcode)
eststo regivi : ivregress 2sls $y (${x} c.${x}#ib3.${m} = ///
${z} c.${z}#ib3.${m}) ib3.$m $c $f , vce(cluster idcode)
* --------------------------------------------------
* 导出回归结果
* --------------------------------------------------
#delimit;
esttab olsbase olsmain olsint regivm regivi using mytable.rtf, append label
title(回归结果) b(3) se(3) brac nogap compress order($x *.${m}*)
note(注：括号中为聚类稳健标准误) coef(_cons "常数项")
indicate("年份=*.year" "行业=*.ind_code",label("控制" "")) nobase interact(" × ")
mgroup("基本模型" "OLS回归" "工具变量回归",pattern(1 1 0 1 0))
mtitle("基本" "主效应" "调节效应" "主效应" "调节效应")
stat(N r2 r2_a F chi2, fmt(0 3 3 3 3) star(F chi2)
label("样本量" "{\i R}²" "调整{\i R}²" "{\i F} 值" "{\i χ}²"));
#delimit cr coef() 选项中成对指定变量的标签，如果不含空格，可以不用双引号，但我个人建议使用双引号来区分变量名和标签。常用用法为 coef(varName1 "varLabel1" varName2 "varLabel2")
结果如下：
带有中文“常数项”的回归结果
这么长一串才进入正文也是独我一份了吧，那么，废话不多说，进入正题，先来看多方程模型的结果导出。
首先多方程模型一般有两种，一种是分步回归（如heckman，2SLS，bdiff），另一种是结构方程模型SEM（如sureg，sem，gsem）。Stata展示结果时涉及两种情况：(1)分步回归一般都是多表格展示结果（或者通常略去第一步的结果）；(2)SEM通常是把多个方程在一张表里展示。对于情况(1)，esttab直接导出一般只能导出最后一步结果，因此需要自己跑一遍前面的步骤；对于情况(2)，esttab可以直接导出结果。
eststo olsbase: reg $y $c $f , vce(cluster idcode)
eststo olsmain: reg $y $x $c $f , vce(cluster idcode)
/* 针对情况(1) 分步回归 Stata通常不汇报第一阶段 */
eststo regivm1: reg $x $z $c $f , vce(cluster idcode) // 自己手动跑第一阶段回归并保存
eststo regivm2: ivregress 2sls $y (${x} = ${z}) $c $f , vce(cluster idcode) first // 指定 first 可以让 Stata 汇报第一阶段结果，但是 esttab 无法直接导出，但可以用来验证是否正确
#delimit ;
esttab olsbase olsmain regivm1 regivm2 using mytable.rtf, replace label
title(回归结果) b(3) se(3) brac nogap compress order($x)
note(注：括号中为聚类稳健标准误) coef(_cons "常数项")
indicate("年份=*.year" "行业=*.ind_code",label("控制" "")) nobase
mgroup("基本模型" "OLS回归" "工具变量回归",pattern(1 1 1 0))
mtitle("基本" "主效应" "第一阶段" "第二阶段")
stat(N r2 r2_a F chi2, fmt(0 3 3 3 3) star(F chi2)
label("样本量" "{\i R}²" "调整{\i R}²" "{\i F} 值" "{\i χ}²"));
#delimit cr结果如下
包含第一阶段的两阶段工具变量回归
下面讲第二种情况：多方程回归，Stata在一张表内展示结果，我们用 heckman 命令来演示（注意heckman的汇报Stata是优先汇报第二阶段，而且其中非常重要的mills比率lambda是在单独的方程中，非常有代表性）
/* 假定我们的 $y 只有 $s == 1 的时候才有取值 */
glo s "south"
// 尽管 heckman 有两种估计，但两种都是在一张表内汇报结果的，在这里一起演示
eststo hecktwo: heckman $y $x $c $f , select($s = $z $c $f ) two
eststo heckmle: heckman $y $x $c $f , select($s = $z $c $f ) mle
/* 1种成熟方案 + 1种折腾但是依旧不好的方案 */
* 成熟方案 手动分步回归存储
* 优点是信息完备 缺点是步骤有些繁琐
est restore hecktwo
eststo heckman1: probit $s $z $c $f if e(sample) // mle 和 two 的第一阶段完全一致不用做两遍了
#delimit;
esttab heckman1 hecktwo heckmle using mytable.rtf, replace label
equation(main = 1, other = .:2:2, IMR = .:3:3)
order(${x}) drop(other:) nonumber noeqline eqlabel("",none)
title(多方程回归结果) b(3) se(3) brac nogap compress
note(注：括号中为聚类稳健标准误) coef(_cons "常数项")
indicate("年份=*.year" "行业=*.ind_code",label("控制" "")) nobase
mtitle("Step 1" "两阶段Heckman" "MLE Heckman")
stat(N chi2, fmt(0 3) star(chi2)
label("样本量" "{\i χ}²"));
#delimit cr
* 穷折腾方案 最后还是不讨好 有大神可以优化
* 优点是不需要跑更多回归了
* 缺点是缺第一阶段统计量；
* 除非愿意先放第二阶段，否则统计量位置不对；
* IMR需要手动标星
#delimit;
esttab hecktwo heckmle using mytable.rtf, replace label
equation(Step2 = 1, Step1 = 2, other = 3) unstack
order(Step1: Step2:${x}) drop(other:) nonumber eqlabel("Step 1" "Step 2")
title(多方程回归结果) b(3) se(3) brac nogap compress
note(注：括号中为聚类稳健标准误) coef(_cons "常数项")
indicate("年份=*.year" "行业=*.ind_code",label("控制" "")) nobase
mtitle("两阶段Heckman" "MLE Heckman")
stat(lambda selambda N chi2, fmt(3 3 0 3) star(chi2)
label("IMR" "S.E. IMR" "样本量" "{\i χ}²"));
#delimit cr
成熟方案结果
主要关注第二和第三行的导出命令：
equation() 用于对每一个多方程模型的方程分类，这里分了三类，main代表主要的，有用的方程，分别是三个模型的第一个方程（main = 1）,other代表要删去的 heckman报告的、缺信息的第一阶段方程，是第二、第三个模型的第二个方程（other = .:2:2）不同模型用英文冒号隔开，如果该模型没有归属这一类的方程则用英文句点表示，IMR是用来衡量Heckman效果的参数（lambda，rho，sigma）在第二、第三个模型的第三个方程（IMR= .:3:3）
drop(other:) 代表扔掉归属 other 的方程内容，也可以用来扔掉变量drop(varName)或者某个模型的变量drop(eqName:varName)
nonumber 由于一个模型只有一个编号，因此不如不要编号
noeqline 将不同方程的分割线去掉
eqlabel("",none) 加上 none 代表不显示方程的标签，若要显示，可参考下文说明
穷折腾方案结果
依然关注第2和第3行命令以及stats()内容
unstack 代表将不同的方程分列显示，因此我们可以在不同列分别显示第一、第二阶段结果
order(Step1: Step2:${x}) 先显示归属 Step1 的方程，在 Step2 方程中优先显示 $x 变量
drop(other:) 扔掉第三个方程内容，用其他方式汇报 IMR
eqlabel("Step 1" "Step 2") 按照方程显示的顺序（可用order指定）分别给方程一个标签
stat(lambda selambda) 用于显示 IMR 的系数和标准误，需要手动据此判断是否需要标星
第二个问题是给变量分组。特别是当我们需要控制很多控制变量时，通过分类能使我们的结果更具有条理。
核心功能就是 refcat() 选项，该选项本身是用于帮助指定离散变量的参考值（比如ib0.foreign 意味着 foreign == 0 为baseline，不参与系数估计）但是活用这个功能可以给我们带来非常漂亮的结果
/* 依然选择之前保存的 olsbase olsmain olsint regivm regivi 回归结果予以说明 */
#delimit;
esttab olsbase olsmain olsint regivm regivi , append label
refcat(grade "{\i\b Controls}" 2.ind_code "{\i\b Fixed Effect}",nolabel)
title(回归结果) b(3) se(3) brac nogap compress order($x *.${m}*)
note(注：括号中为聚类稳健标准误) coef(_cons "常数项")
indicate("年份=*.year" ,label("控制" "")) nobase interact(" × ")
mgroup("基本模型" "OLS回归" "工具变量回归",pattern(1 1 0 1 0))
mtitle("基本" "主效应" "调节效应" "主效应" "调节效应")
stat(N r2 r2_a F chi2, fmt(0 3 3 3 3) star(F chi2)
label("样本量" "{\i R}²" "调整{\i R}²" "{\i F} 值" "{\i χ}²"));
#delimit cr
带上分组的控制变量 表格过长没有放全
关键看esttab第二行的 refcat() 用法
refcat(grade "{\i\b Controls}" 2.ind_code "{\i\b Fixed Effect}",nolabel) 表示在变量 grade 和 2.ind_code 上方加一行 加粗斜体 的 和 。2.ind_code 是因为默认把 1.ind_code作为行业的参考值。注意一定要加上 ,nolabel 这个 suboption
终于到了最后的 矩阵导出 了
参考连玉君老师 [一文读懂中介效应分析](https://zhuanlan.zhihu.com/p/99435552)的内容，使用bootstrap和sgmediation来进行检验时，最终的结果是通过 estat boot 导出的，并不能直接输出到word文档中。
bootstrap r(ind_eff) r(dir_eff) r(tot_eff), reps(500) : ///
sgmediation ${y}, mv(${m}) iv(${x}) cv(${c} ${f}) //计算中介效应和间接效应
estat bootstrap, bc //计算置信区间 由于 percentile 参考价值不高故只报告 bc 的置信区间
/* estat boot, bc 结果如下 */
Bootstrap results Number of obs = 28,104
Replications = 500
Command: sgmediation2 ln_wage, mv(race) iv(hours) cv(grade age i.year i.ind_code)
_bs_1: r(ind_eff)
_bs_2: r(dir_eff)
_bs_3: r(tot_eff)
------------------------------------------------------------------------------
| Observed Bootstrap
| coefficient Bias std. err. [95% conf. interval]
-------------+----------------------------------------------------------------
_bs_1 | -.00019584 -1.40e-07 .00002278 -.0002436 -.0001561 (P)
| -.0002438 -.0001564 (BC)
_bs_2 | .00103865 .0000167 .0003553 .0004068 .0017334 (P)
| .0004121 .001745 (BC)
_bs_3 | .00084281 .0000166 .00035402 .0002039 .0015258 (P)
| .0001944 .001519 (BC)
------------------------------------------------------------------------------
Key: P: Percentile
BC: Bias-corrected我尝试对结果进行调整，最终结果如下
间接效应路径系数与置信区间
代码如下
bootstrap r(ind_eff) , reps(500) : /// 由于中介效应是否显著仅关注简介效应就行，为此我们只汇报 r(ind_eff)
sgmediation2 ${y}, mv(${m}) iv(${x}) cv(${c} ${f})
estat boot, bc percentile // 查看结果方便对照
cap mat drop Boot
mat Boot = [e(b)',e(ci_bc)']
local lx: var label $x
local lmed: var label $m
local ly: var label $y
mat rownames Boot = "`lx' -> `lmed' -> `ly'" // 多中介变量以空格分隔，依次命名
mat colnames Boot = "路径系数" "95% 置信区间" // 最后两列在 word 中合并单元格即可
esttab mat(Boot, fmt(a3 a3 a3)) using mytable2.rtf, replace nomtitle核心在于 用 esttab 导出矩阵
系数和置信区间分别存储在 e(b) 和 e(ci_bc) 中，加英文单引号'是为了赚置矩阵
mat Boot = [e(b)',e(ci_bc)'] 将路径系数与置信区间放在一个矩阵中
mat rownames 和 mat colnames 分别指定了矩阵的名称
好啦，真的被榨干了，一滴都不剩了。esttab其实还有很多可以玩的花活，不过不遇到还真想不起来还有啥冷门的导出。
希望写的足够清晰能让大家看明白
