Stata：手动计算置信区间

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

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

🍓 课程推荐：2024 机器学习与因果推断专题
主讲老师：司继春 (上海对外经贸大学) ；张宏亮（浙江大学）
课程时间：2024 年 11 月 9-10 日 ；16-17日
课程咨询：王老师 18903405450（微信）

课程特色 · 2024机器学习与因果推断：

• 懂原理、会应用。本次课程邀请了两位老师合作讲授，目的在于最大限度地实现理论与应用的有机结合。为期四天的课程，分成两个部分：第一部分讲解常用的机器学习算法和适用条件，以及文本分析和大语言模型；第二部分通过精讲 4-6 篇发表于 Top 期刊的论文，帮助大家理解各类机器学习算法的应用场景，以及它们与传统因果推断方法的巧妙结合。
• 以 Top 期刊论文为范例。目前多数人的困惑是不清楚如何将传统因果推断方法与机器学习结合起来。事实上，即便是 MIT 和 Harvard 的大牛们也都在「摸着石头过河」。为此，通过论文精讲和复现来学习这部分内容或许是目前最有效的方式了。张宏亮老师此前在浙江大学按照这一模式教授了「因果推断和机器学习」课程，效果甚佳：学生们能够逐渐建立起研究设计的理念，并在构造识别策略时适当地嵌入机器学习方法。

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

作者：刘东 (中国农业大学)
邮箱：ld01@cau.edu.cn

编者按：本文主要参考自以下文章，特此致谢！

• Mark Bounthavong. 2022. Blog. A tutorial on how to add the 95% CI to a two-way line plot. -Link-
• Jann B. Plotting regression coefficients and other estimates[J]. The Stata Journal, 2014, 14(4): 708-737. -PDF-

1. 均值的置信区间：手动操作

1.1 折线图

在描述性统计中，常常需要展示某些变量的年度均值及其时序特征，有时还需要附加显示均值的 95% 置信区间。在这个示例中，我们将使用 Stata 手册中的数据文件 nlswork.dta 进行演示。首先，我们将调入数据，并生成支出变量的年度均值，然后展示这些均值随年份变化的趋势：

. webuse "nlswork.dta", clear
. replace year = 1900 + year
. tab year 
. gen wage = exp(ln_wage)
. label var wage "hourly wage/GNP deflator"
. egen wage_yr = mean(wage), by(year)
. sort year
. graph twoway connected wage_yr year,  ///
>     title("Average hourly wage of US young women (1968-1988)") ///
>     ytitle("Average hourly wage ($)") ///
>     xtitle("Year")                    ///
>     xlab(1968(1)1988, angle(30))      ///
>     ylab(, labsize(small) nogrid)     ///
>     graphregion(color(white))         ///
>     bgcolor(white)
. graph export "nlswork_connected1.png", width(600) replace

输出结果如下：

1.2 添加置信区间

上图仅呈现了各年度的平均小时工资水平，但无法知晓各个年度上的平均工资的标准误。我们可以使用 ci mean 命令来估计年度均值的 95% 置信区间 (CI)。

. webuse "nlswork.dta", clear
. replace year = 1900 + year
. tab year 
. gen wage = exp(ln_wage)
. label var wage "hourly wage/GNP deflator"
. egen wage_yr = mean(wage), by(year)

. * 添加置信区间 
. gen high = .
. gen low = .
. * 为每年的均值生成置信区间
. forvalues i = 1968/1988 {
  2.         ci mean wage if year == `i'
  3.         replace high = r(ub) if year == `i'
  4.         replace  low = r(lb) if year == `i'
  5. }
. sort year
. graph twoway                       ///
>         (rcap low high year)       ///
>         (connected wage_yr year,   ///
>         color(navy) msize(small)   ///
>         title("Average hourly wage of US young women (1968-1988)") ///
>         ytitle("Average hourly wage ($)")           ///
>         xtitle("Year")                              ///
>         xlab(1968(1)1988, angle(30) labsize(small)) ///
>         ylab(, labsize(small) nogrid) ///
>         graphregion(color(white))     ///
>         bgcolor(white)    ///
>         color(navy)       ///
>         legend(ring(0) position(5)))
. graph export "nlswork_connected_CI_bars.png", width(600) replace

上述代码中最核心的是如下三行：

ci mean wage if year == `i'
replace high = r(ub) if year == `i'
replace low = r(lb) if year == `i'

这段代码的目的是计算每一年 (1968 至 1988 年) 的工资均值的置信区间，并将置信区间的上下限分别存储在变量 high 和 low 中。具体说明如下：

计算置信区间：

ci mean wage if year == `i'

这行代码使用 ci mean 命令计算特定年份 i 的小时工资 wage 的均值及其置信区间。ci mean 命令根据数据生成置信区间，并将结果存储在返回值 r(ub) (上置信限) 和 r(lb) (下置信限) 中。置信区间的计算公式为：

其中： 是样本均值； 是 分布的临界值 (自由度为 ； 是样本标准差； 是样本量)。

存储置信区间的上下限：

replace high = r(ub) if year == `i'
replace low = r(lb) if year == `i'

这两行代码分别将计算得到的置信区间上限值 r(ub) 和下限值 r(lb) 存储到变量 high 和 low 中，条件是当前年份等于循环变量 i。

通过这些操作，我们可以得到每一年工资均值的置信区间，为后续的分析和图形展示提供更详细的信息。置信区间有助于评估均值的可靠性和波动范围，使分析结果更具说服力。

2. 更快捷地呈现均值和置信区间：coefplot 命令

采用上文介绍的 ci mean 命令虽然也可以达到计算并呈现置信区间的目的，但需要编写较多的代码。尤其是当我们需要进行多组样本的对比分析时，上述方法就会变得很繁琐。

此时，我们可以使用 Bann 编写的 coefplot 命令。该命令经常出现在 AER，QJE，JPE，JFE 等 Top 期刊论文的复现文档中。我们可以使用 ssc install 命令安装 coefplot 命令：

. ssc install coefplot, replace 
. help coefplot

相比于 Stata 自带的 [R] marginsplot 命令，coefplot 命令更加灵活，可以实现更为复杂的作图需求。

2.1 coefplot：单模型系数可视化

. sysuse auto, clear
. regress price mpg trunk length turn
. coefplot, drop(_cons) xline(0)

需要说明的是，coefplot 不仅可以在执行 regress 命令后对估计系数进行可视化呈现，它同样适用于 logit，probit 等命令。

2.2 coefplot：多模型系数可视化对比

在实证分析中，常常需要将多个模型的估计结果进行对比。此时，我们可以借助 coefplot 命令来绘制多模型的回归系数图。下面的例子展示了如何使用 coefplot 命令来对比国内汽车和国外汽车的价格回归模型。

2.2.1 合图呈现两个模型的系数

. sysuse auto, clear
. regress price mpg trunk length turn if foreign==0
. estimates store D
. regress price mpg trunk length turn if foreign==1
. estimates store F
. * 要为单个图指定单独的选项，将该模型名称及其选项扩在括号中
. coefplot (D, label(Domestic Cars) pstyle(p3)) ///
>          (F, label(Foreign Cars) pstyle(p4)), ///
>          msymbol(S) drop(_cons) xline(0)

输出结果如下：

2.2.2 两个子图

我们也可以把上图拆分成两个子图来呈现：

. coefplot D, bylabel(Domestic Cars) ||  ///
>          F, bylabel(Foreign Cars)  ||, ///
>          drop(_cons) xline(0)

2.2.3 多模型子图

. * 国产车
. regress weight mpg trunk length turn if foreign == 0
. estimates store wD

. * 进口车
. regress weight mpg trunk length turn if foreign == 1
. estimates store wF

. * 使用 coefplot 绘制回归系数图
. coefplot (D, label(Domestic)) (F, label(Foreign)), ///
>     bylabel(Price) || wD wF, ///
>     bylabel(Weight) ||,      ///
>     drop(_cons) xline(0)     ///
>     byopts(xrescale)

2.2.4 比较单变量回归和多元回归系数差异

. sysuse auto, clear
. * 多元回归模型
. regress price mpg trunk length turn
. estimates store multivariate

. * 循环对每个变量进行单变量回归
. foreach var in mpg trunk length turn {
  2.     quietly regress price `var'
  3.     estimates store `var'
  4. }

. * 使用 coefplot 绘制回归系数图
. coefplot (mpg \ trunk \ length \ turn, label(bivariate)) ///
>          (multivariate), ///
>          drop(_cons) xline(0)

2.2.5 匹配系数和方程

coefplot 的默认设置是使用每个模型中的第一个 (非零) 方程，并通过它们的名称 (忽略方程名) 来匹配系数。例如，regress 返回一个 (无名称的) 方程，其中包含回归系数，而 tobit 返回两个方程，一个名为 model 的方程包含回归系数，另一个名为 sigma 的方程包含回归的标准误差。因此，coefplot 的默认设置是匹配两个模型中的回归系数，并忽略 tobit 模型中的方程 sigma。

. sysuse auto, clear
. foreach v of var price mpg trunk length turn {
  2.         quietly sum `v'
  3.         quietly replace `v' =(`v' - r(mean)) / r(sd)
  4. }
. regress price mpg trunk length turn
. estimate store regress
. tobit price mpg trunk length turn, ll(-.5)
. estimate store tobit
. coefplot regress tobit, xline(0)

2.3 设置多个置信区间

coefplot 的默认设置是绘制 95% 置信区间。要指定不同的置信水平或包含多个置信区间，可以使用 levels() 选项。下面是一个包含 99.9%、99% 和 95% 置信区间的示例：

. sysuse auto, clear
. regress price mpg trunk length turn
. coefplot, drop(_cons) xline(0) msymbol(s) mfcolor(white) levels(99.9 99 95) ///
>     legend(order(1 "99.9" 2 "99" 3 "95") row(1))

如果估计命令提供了预先计算的置信区间，可以使用 ci() 选项将它们包含在绘图中。例如，要绘制通过 bootstrap 在 e(ci normal)、e(ci percentile) 和 e(ci bc) 中提供的正态近似、百分位和偏差校正置信区间，可以输入以下命令：

. regress price mpg trunk length turn, vce(bootstrap)
. coefplot (, ci(ci_normal) label(normal)) (, ci(ci_percentile) label(percentile)) ///
>          (, ci(ci_bc) label(bc)), drop(_cons) xline(0) legend(row(1))

除了 levels() 和 ci() 之外，您还可以使用选项 cismooth 添加平滑的置信区间。默认情况下，cismooth 为 50 个等间距水平 (1、3、...、99) 生成置信区间，并使用渐变的颜色强度和不同的线宽，如下面的示例所示：

. coefplot, drop(_cons) xline(0) cismooth grid(none)

2.4 其他的绘图风格

垂直模式，默认情况下，coefplot 生成一个水平图形，y 轴上有标签，x 轴上有值。要翻转轴，请指定垂直选项。

. sysuse auto, clear
. regress price mpg trunk length turn
. coefplot, drop(_cons) vertical yline(0)

要更改用于标记和置信区间的绘图类型，可以使用 recast() 选项。可用于标记的绘图类型包括标准的 twoway 图，例如 scatter (默认)、line、dot 或 bar。对于置信区间，可以使用范围图，例如 rspike (默认)、rline、rcap 或 rbar。

. coefplot, drop(_cons) xline(0) ciopts(recast(rcap))

要将系数的值作为标记标签添加到图表中，可以使用 mlabel 选项，还可以结合使用 format() 来设置显示格式。

. sysuse auto, clear
. keep if rep78>=3
. regress mpg headroom i.rep##i.foreign
. coefplot, xline(0) mlabel format(%9.2g) mlabposition(12) mlabgap(*2)

另一种显示数值标签的方法：

. sysuse auto, clear
. keep if rep78 >= 3
. regress mpg headroom i.rep##i.foreign
. mata: st_matrix("e(box)", (st_matrix("e(b)") :- 2 \ st_matrix("e(b)") :+ 2))
. coefplot,    ///
>     xline(0) ///
>     mlabel format(%9.2g) ///
>     mlabposition(0)      ///
>     msymbol(i)           ///
>     ci(95 box)           ///
>     ciopts(recast(. rbar) barwidth(. 0.35) fcolor(. white) lwidth(. medium))

3. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 置信区间, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：其它
• 劳伟健, 2022, Stata：正确理解置信区间, 连享会 No.1086.
• 专题：Stata程序
• 吴雄, 童天天, 2020, Stata：Bootstrap-自抽样-自举法, 连享会 No.43.
• 专题：回归分析
• 曹琳君, 2021, Stata：如何估计置信区间？, 连享会 No.712.
• 专题：Stata入门
• 肖蕊, 2020, 25常见种误区：P值、置信区间和统计功效, 连享会 No.221.
• 专题：Stata绘图
• 雷诺, 2023, Stata绘图：confcomptwp-二维系数的置信区间和可比区间可视化, 连享会 No.1257.
• 马洪栋, 2024, Stata绘图：高级柱状图（一）-均值和置信区间-cibar-coefpl, 连享会 No.1379.

🍓 课程推荐：2024 机器学习与因果推断专题
主讲老师：司继春 (上海对外经贸大学) ；张宏亮（浙江大学）
课程时间：2024 年 11 月 9-10 日 ；16-17日
课程咨询：王老师 18903405450（微信）

尊敬的老师 / 亲爱的同学们：

连享会致力于不断优化和丰富课程内容，以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求，我们诚挚地邀请您参与到我们的课程规划中来。请您在下面的问卷中，分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源，将直接影响到我们课程的改进和创新。我们期待您的反馈，因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间，与我们共同塑造更加精彩的学习旅程！https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见！

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

• 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。