Stata绘图：COVID-19数据可视化-山脊图

👇 连享会 · 推文导航 | 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 的大牛们也都在「摸着石头过河」。为此，通过论文精讲和复现来学习这部分内容或许是目前最有效的方式了。张宏亮老师此前在浙江大学按照这一模式教授了「因果推断和机器学习」课程，效果甚佳：学生们能够逐渐建立起研究设计的理念，并在构造识别策略时适当地嵌入机器学习方法。

作者：王胜文 (山东财经大学)
邮箱：sw8258@foxmail.com

编者按：本文整理自 COVID-19 visualizations with Stata Part 8: Ridgeline plots (Joy plots) ，作者 Asjad Naqvi，特此致谢！

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

目录

• 1. 前言

• 2. 任务一：数据设置

• 3. 任务二：按顺序获取面积图

• 4. 任务三：图形的全自动化

• 5. 牛刀小试

• 6. 其他的 Stata 指南

• 7. 作者 Asjad Naqvi 的简介

• 8. 相关推文

在本操作指南中，我们将要学习如何从 Our World in Data 网站中获取公开的新冠肺炎数据，并在 Stata 中绘制“山脊图” (也叫“峰峦图”或“堆叠图”)。在本指南的最后，我们将学习绘制下面的图片。

1. 前言

本指南是假设您已具备 Stata 操作技能的基础上所撰写的。如果您是 Stata 新手，或者第一次使用本指南，那么强烈建议您先阅览 指南1 与 指南2。和之前几个操作指南一样，标准化介绍也适用于本指南。我们先从工作流程管理的文件夹结构开始介绍。

在名为 graphs 的文件夹中，我们创建了一个名为 guide8 的子文件夹，用于存储操作过程中生成的图形。有关如何构建您文件夹的详细信息，请参考 指南1。为了使图片与此处所示完全一致，还需要设置几个附加命令：

• Cleanplots主题 (指南2)

net install cleanplots, from("https://tdmize.github.io/data/cleanplots")
set scheme cleanplots, perm

• 调色板套装 (指南2)

net install palettes, replace from("https://raw.githubusercontent.com/benjann/palettes/master/")
net install colrspace, replace from("https://raw.githubusercontent.com/benjann/colrspace/master/")

• 将默认图形字体设置为 Arial Narrow

graph set window fontface "Arial Narrow"

另一种方法是生成一些图形，您可以右键单击图形窗口，单击属性，然后选择图形的字体。这一步是可选的，我将在另一个指南中介绍 Stata 中字体的使用。

2. 任务一：数据设置

在指南1与指南2中，讨论了从互联网中提取数据的详细过程。但是为了本指南的完整性，此处提供了数据获取的基本代码。首先，我们在 Our World in Data COVID-19 repository 中获取数据并设置数据。

*** COVID 19 data ***
insheet using "https://covid.ourworldindata.org/data/owid-covid-data.csv", clear
save ./raw/full_data_raw.dta, replace
gen date2 = date(date, "YMD")
format date2 %tdDD-Mon-yy
drop date
ren date2 date
ren location country
replace country = "Slovak Republic" if country == "Slovakia"
drop if date < 21915  // 1st Januarycompress
save "./master/OWID_data.dta", replace

接下来，我们需要保留一个国家样本，并稍微清理一下数据。请注意，此处可以使用任何国家的数据集。

use ./master/OWID_data.dta, clear
summ date
drop if date >= `r(max)' - 2       // drop the last two observations
gen group = .
replace group = 1 if               ///
    country == "Austria"     |     /// 
    country == "Belgium"     |     ///
    country == "Czech Republic"  | ///
    country == "Denmark"     |     ///
    country == "Finland"     |     ///
    country == "France"      |     ///
    country == "Germany"     |     ///
    country == "Greece"      |     ///
    country == "Hungary"     |     ///
    country == "Italy"       |     ///
    country == "Ireland"     |     ///
    country == "Netherlands" |     ///
    country == "Norway"      |     ///
    country == "Poland"      |     ///
    country == "Portugal"    |     ///
    country == "Slovenia"    |     ///
    country == "Slovak Republic" | ///
    country == "Spain"       |     ///
    country == "Sweden"      |     ///
    country == "Switzerland" |     ///
    country == "United Kingdom"
keep if group==1
keep country date new_cases new_deaths
 // generate a numerical country var
encode country, gen(country2)  
order country country2 date

*** fix some data errors
replace new_cases  = 0 if new_cases < 0
replace new_deaths = 0 if new_deaths < 0
lab var new_cases "New cases"
lab var new_deaths "New deaths"
lab var date "Date"
drop if date < 21960  
// drop data points before 1st Jan 2020
format date %tdDD-Mon
des

最后我们应该只剩下五个变量，如下所示：

ontains data from ./master/OWID_data.dta
 Observations:        19,742                  
    Variables:             5                  8 Nov 2022 14:44
---------------------------------------------------------------
Variable      Storage   Display    Value
    name         type    format    label      Variable label
---------------------------------------------------------------
country         str32   %32s                  
country2        long    %15.0g     country2   
date            float   %tdDD-Mon             Date
new_cases       long    %12.0g                New cases
new_deaths      int     %8.0g                 New deaths
---------------------------------------------------------------
Sorted by: 
     Note: Dataset has changed since last saved.

3. 任务二：按顺序获取面积图

在这部分，我们没有像在上篇指南那样使用 7 天滑动平均值，而是使用了非参数 Lowess 平滑函数。举个例子，country2==1 的 Lowess 拟合可以通过下列代码生成：

twoway (lowess new_cases date if country2==1, bwid(0.1)) ///
  (scatter new_cases date if country2==1, msize(small) mcolor(gs8%30)), legend(off)
graph export ./graphs/guide8/graph0.png, replace wid(5000)

请注意，我只是将图形导出命令作为例子在这里展示，如果您没有使用上述文件夹结构，那么此行代码则起不到作用，可以删除。通过上面的代码我们可以得到下面的图形。

在上面的代码中，bwid(0.1) 是带宽为 0.1 的缩写，并使用默认的 Epanechnikov 核函数。数值越小，拟合越接近实际数据；而当数值越大，拟合线就越平滑。对于非参数函数来说没有“绝对合适的”带宽，所以就需要具体问题具体分析。我们先从前四个国家被定义的 country2 开始绘制：

twoway (lowess new_cases date if country2==1, bwid(0.1)) ///
    (lowess new_cases date if country2==2, bwid(0.1))    ///
    (lowess new_cases date if country2==3, bwid(0.1))    ///
    (lowess new_cases date if country2==4, bwid(0.1)), legend(off)
graph export ./graphs/guide8/graph1.png, replace wid(2000)

该图以平滑线图的形式显示了每日报告的病例：

对于山脊图来说，变量一般是经过标准化的 (通常在 0 到 1 之间)，要么是每个类别的最大值 (在我们的例子中是国家)，要么是根据感兴趣变量数据中的总体最大值。前者更容易比较一个国家随时间变化发生的变化，而后者更适用于比较所有国家之间产生的变化。

如果差异很大，那么“全球”标准化可能会隐藏规模较小的国家。如上图所示，黑线是蓝线很小的一部分。我们可以使用每个与人口相关的案例来更准确地比较各个国家和地区的不同，但是对于本指南，我们坚持以“地区”一级进行国家层面的标准化。上述过程可通过以下方式实现。

***每个国家的病例标准化过程 (范围为0-1) 
gen cases_norm = . 
gen deaths_norm = . 
levelsof country2, local(levels)
foreach x of local levels{
    summ new_cases if country2==`x'
    replace cases_norm = new_cases / `r(max)' if country2==`x'
    summ new_deaths if country2==`x'
    replace deaths_norm = new_deaths / `r(max)' if country2==`x'
 }

请注意使用levelsof和locals生成新的变量。一旦有了标准化变量，我们就可以将较小的拟合值存储为单独变量，这有助于我们将数据导进随后的图形中：

lowess cases_norm date if country2==1, bwid(0.1) gen(ytop1) nograph
lowess cases_norm date if country2==2, bwid(0.1) gen(ytop2) nograph
lowess cases_norm date if country2==3, bwid(0.1) gen(ytop3) nograph
lowess cases_norm date if country2==4, bwid(0.1) gen(ytop4) nograph

对于每个国家，我们生成一个新变量 ytopx ，其中包含最低平滑拟合值。我们还可通过以下命令绘制图片。

twoway (line ytop1 date) (line ytop2 date) (line ytop3 date) ///
    (line ytop4 date), legend(off)

我们可以绘出下面的图，其中 Y 轴上数值的范围是 0 到 1 。为了适应此数据范围，所有国家的曲线都已按比例进行了调整。从技术上看，如果平滑函数的带宽 (目前设置为 0.1 ) 减少，那么线条也将在 Y 轴的上限上接近 1，因为相邻观测值的权重较小。此图能能加清晰地显示出每个国家之间的内部差异：

我们也可以将这些线条转化为面积，在这里我们使用 twoway area 函数，这个函数需要输入三个变量，分别是一个 X 轴值和两个 Y 轴的值分别作为上限和下限。通过下列代码可以完成此步骤：

gen y0 = 0
twoway (rarea ytop1 y0 date) (rarea ytop2 y0 date) ///
    (rarea ytop3 y0 date) (rarea ytop4 y0 date), legend(off)

可以生成下图：

为使图片更加美观，我们可以使用下面代码调整填充色彩 (需使用 Stata.15 及更高版本)。

twoway (rarea ytop1 y0 date, fcolor(%60)) ///
  (rarea ytop2 y0 date, fcolor(%60))      ///
  (rarea ytop3 y0 date, fcolor(%60))      ///
  (rarea ytop4 y0 date, fcolor(%60)), legend(off)

通过上述代码，我们可以得到一个显示所有国家趋势的面积图：

如果我们通过下面代码颠倒各个国家的顺序：

twoway (rarea ytop4 y0 date, fcolor(%80)) ///
  (rarea ytop3 y0 date, fcolor(%80))      ///
  (rarea ytop2 y0 date, fcolor(%80))      ///
  (rarea ytop1 y0 date, fcolor(%80)), legend(off)

我们可以得到以下图片：

可以看出，绘制图片的顺序被颠倒了，这一步对我们创建山脊图来说非常重要，因为需要控制绘图顺序 (先看规则、最后绘制图片)。同时，这一步对控制颜色方案也很重要，我们将稍后进行讨论。

在下一步中，我们将向图中添加两个控制因素。由于每个国家都在 0 到 1 之间进行了标准化，因此我们只生成底部的 Y 值，作为从 1 到国家数量的数字 (本例中国家数量为 4 ) 。同时我们也移动了顶部的部分，这些部分显示了数值的分布：

gen y1 = 1
gen y2 = 2
gen y3 = 3
gen y4 = 4
gen ytest4 = ytop4 + y4
gen ytest3 = ytop3 + y3
gen ytest2 = ytop2 + y2
gen ytest1 = ytop1 + y1
twoway (rarea ytest4 y4 date, fcolor(%80)) ///
  (rarea ytest3 y3 date, fcolor(%80))      ///
  (rarea ytest2 y2 date, fcolor(%80))      ///
  (rarea ytest1 y1 date, fcolor(%80)), legend(off)

通过上面的代码我们可得到下图：

所有国家都堆叠在一起，并且区间都在 0 到 1 之间，这基本上是生成山脊图所需要的核心图形。现在我们要将 Y 轴压扁，只需要在面积图的“底端”部分生成一个变量，从而将 Y 轴刻度之间的空间减少四分之一。

**** squish the axis 
gen y1_new = 1 / 4 
gen y2_new = 2 / 4
gen y3_new = 3 / 4
gen y4_new = 4 / 4
twoway (rarea ytest4 y4_new date, fcolor(%80)) ///
  (rarea ytest3 y3_new date, fcolor(%80))      ///
  (rarea ytest2 y2_new date, fcolor(%80))      ///
  (rarea ytest1 y1_new date, fcolor(%80)), legend(off)

如果我们将绘制的图片使用原始的“顶部”层，我们会得到如下图片：

同时，我们还需要向下移动顶部分散层，具体操作如下：

**** and shift the top curves  
gen ytest1_new = ytop1 + y1_new
gen ytest2_new = ytop2 + y2_new
gen ytest3_new = ytop3 + y3_new
gen ytest4_new = ytop4 + y4_new
twoway (rarea ytest4_new y4_new date, fcolor(%80)) ///
  (rarea ytest3_new y3_new date, fcolor(%80)) ///
  (rarea ytest2_new y2_new date, fcolor(%80)) ///
  (rarea ytest1_new y1_new date, fcolor(%80)), legend(off) 
graph export ./graphs/guide8/graph8.png, replace wid(2000)

由此我们可以得到这个图：

要注意的是，Y 轴压缩了 4 倍，我们也可以通过调整倍数来获取最优的结果。下一步我们生成用于标记国家的散点。

// labels for the country namesegen tag = tag(country)
summ date
gen xpoint = `r(min)' if tag==1gen ypoint = .
replace ypoint = (ytest1_new + 1/20) if xpoint!=. & country2==1 
replace ypoint = (ytest2_new + 1/20) if xpoint!=. & country2==2 
replace ypoint = (ytest3_new + 1/20) if xpoint!=. & country2==3 
replace ypoint = (ytest4_new + 1/20) if xpoint!=. & country2==4

这一步很重要，因为国家不在由 Y 轴上的分类变量 (1、2、3、4...) 而定义而是通过 0.25、0.5、0.75 这一类的数进行定义。这些步骤也会根据我们对 Y 轴的挤压程度不同而发生变化。

我们没有试图找出 Y 轴上需要标记的值，而是完全去掉了 Y 轴，并用散点图进行了标记。X 轴的值可以简单看作为早期的日期，我们将其标记为 xpoint；而生成面积图的 Y 轴，我们将其标记为 ypoint。同时我们要将 ypoint 向上微调，以避免线条之间的重叠。目前，我们已经生成了标记图片的所有元素：

summ date 
local x1 = `r(min)'
local x2 = `r(max)'
twoway (rarea ytest4_new y4_new date, fcolor(%80)) ///
  (rarea ytest3_new y3_new date, fcolor(%80))      ///
  (rarea ytest2_new y2_new date, fcolor(%80))      ///
  (rarea ytest1_new y1_new date, fcolor(%80))      ///
  (scatter ypoint xpoint, mcolor(white) msize(zero) msymbol(point) ///
  mlabel(country) mlabsize(*0.6) mlabcolor(black)),                ///
  xlabel(`x1'(10)`x2', nogrid labsize(vsmall) angle(vertical))     ///
  ylabel(, nolabels noticks nogrid) yscale(noline) ytitle("") xtitle("") ///
  legend(off)
graph export ./graphs/guide8/graph9.png, replace wid(2000)

同时，我们还需要进行以下工作。首先，我们需要通过局部变量自动地控制数据范围；接下来，我们要添加一个散点，使用隐藏的标记 msize(zero) 做为国家标签。同时我们也完全隐藏了 Y 轴，图中所有的刻度线也都关闭。使用上述代码，我们可以得到下面的图形：

这为我们能够生成自定义颜色与格式的图片提供了便利。

4. 任务三：图形的全自动化

在这一部分中，我们将使用以下几个操作步骤，这些操作步骤的逻辑已经在以前的指南中探讨过。由于我们使用了几个局部变量代码，它们必须要一次性执行，所以我们在这里给出完整代码块，需要整体运行。

cap drop y*   // drop the variables we generated earlier
gen ypoint=.
levelsof country2, local(levels)
local items = `r(r)' + 6
foreach x of local levels {
    summ country2
    local newx = `r(max)' + 1 - `x'   // reverse the sorting
    lowess cases_norm date if country2==`newx', bwid(0.05) gen(y`newx') nograph
    gen ybot`newx' =  `newx'/ 4       // squish the axis
    gen ytop`newx' = y`newx' + ybot`newx'
    colorpalette matplotlib autumn, n(`items') nograph
    local mygraph `mygraph' rarea  ytop`newx' ybot`newx' date, ///
        fc("`r(p`newx')'%75") lc(white) lw(thin) 
    replace ypoint= (ytop`newx' + 1/16) if xpoint!=. & country2==`newx' 
}
summ date 
local x1 = `r(min)'
local x2 = `r(max)'
twoway  `mygraph' (scatter ypoint xpoint, mcolor(white) msize(zero)  ///
    msymbol(point) mlabel(country) mlabsize(*0.6) mlabcolor(black)), ///
    xlabel(`x1'(10)`x2', nogrid labsize(vsmall) angle(vertical))     ///
    ylabel(, nolabels noticks nogrid) yscale(noline)  ytitle("")     ///
    xtitle("") legend(off)                                           ///
    title("{fontface Arial Bold:COVID-19 daily cases in Europe}")    ///
    note("Data sources: Our World in Data, JHU, ECDC.", size(tiny)) 
graph export ./graphs/guide8/graph10.png, replace wid(2000)

现在让我们来仔细分析上面的代码。levelsof country2, local(levels) 将所有国家的特征信息保存在以 levels 未命名的存储块中。通过这个命令，我们还可以弄清楚 levels 的具体数值为多少，具体公式为 local items = r(r) + 6。

同时，我们用这条命令 colorpalette matplotlib autumn, n(items) nograph 来生成新的颜色，在这里为了让色彩更加明显，我还特意添加了额外的 6 个值来扩展颜色范围。

summ country2 和 local newx = r(max)+1-x 这两个命令重置了国家的顺序，使其从上到下依次绘制。虽然有其他方法也可以实现这步操作，但是这两个命令非常有效。lowess 将每个国家级的平滑值存储在国家特定变量中。当然这部分代码也可以通过使用临时变量 tempvar 进行优化，但我们这么做的目的就是为了弄清楚每个步骤生成的数据是什么样的。

在下一步操作中，我们生成顶部与底部的 Y 值，以便生成新的区域。在 gen ybot newx = newx / 4 代码中，数值 4 是“挤压”因子。人们可以根据国家的数量、波动的期望值以及它们之间的期望差距调整这个数字，以获得图片的最佳效果。分母的数字越高，图中的山脊的幅度就越大。

变量 ypoint 生成国家标签的散布。由于我们关闭了坐标轴，所以我们只能使用 X-Y 值来标记每个山脊。我们添加 +1/16 的部分只是将散射点稍微向上移动，以便它们与直线对齐，使图片达到最佳效果。否则，标签将在每条线的中点结束。

在标题中 fontface Arial Bold:xxx 允许我们自定义标题的字体 (有关更多信息，请参阅 Stata 字体指南)。

这幅图清晰显示了欧洲样本国家第二波疫情浪潮的形成。

5. 牛刀小试

尝试生成每日病亡的统计图片，图片使用 colorpalette matplotlib winter 配色方案来完成。

还可以尝试生成每个人的病例和病亡图片，并使用全局标准化。

6. 其他的 Stata 指南

• 第一部分：数据设置和自定义图形简介
• 第二部分：自定义颜色方案
• 第三部分：热力图
• 第四部分：地图
• 第五部分：堆积面积图
• 第六部分：动画图
• 第七部分：双倍时间图
• 第八部分：山脊图
• 第九部分：自定义条形图
• 第十部分：流图

如果您喜欢上述指南并发现它们确实好用，请多多关注 Stata 指南。另外，如果您使用这些指南，请尽量分享、推广您的可视化成果。

7. 作者 Asjad Naqvi 的简介

作者 Asjad Naqvi 是一名专业经济学家，自 2003 年以来一直在使用 Stata。Asjad Naqvi 目前居住在奥地利维也纳，在维也纳经济商业大学 (WU) 和国际应用系统分析研究所 (IIASA) 工作。您可以在 ResearchGate 和 Google Scholar 上找到 Asjad Naqvi 的研究成果，在 GitHub 上找到 Stata 代码库。

Asjad Naqvi 的博客中有关 Stata 内容的链接：Stata 指南，这里会定期发布新的精彩内容。

8. 相关推文

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

• 专题：Stata绘图
• Stata绘图：addplot-层层叠加轻松绘图
• Stata 绘图：用 Stata 绘制一打精美图片-schemes
• 常用科研统计绘图工具介绍
• Stata空间计量：莫兰指数绘图moranplot命令介绍
• Stata绘图-组间差异可视化：不良事件火山图、点阵图
• Stata绘图极简新模板：plotplain和plottig-T251
• 给你的图形化个妆：Stata绘图常用选项汇总-上篇
• 给你的图形化个妆：Stata绘图常用选项汇总-下篇
• Stata绘图：柱状图专题-T212
• Stata绘图：回归系数可视化-论文更出彩
• Stata绘图：世行可视化案例-条形图-密度函数图-地图-断点回归图-散点图
• Stata绘图：随机推断中的系数可视化
• Stata绘图：重新定义坐标轴刻度标签
• Stata绘图：用-bytwoway-实现快速分组绘图
• Stata绘图：一个干净整洁的-Stata-图形模板qlean
• Stata绘图：怎么在Stata图形中附加水平线或竖直线？
• Stata绘图：在图片中添加虚线网格线
• Stata绘图：制作教学演示动态图-GIF
• Stata绘图：绘制一颗红心-姑娘的生日礼物
• Stata绘图：让图片透明——你不要掩盖我的光芒
• Stata绘图：bgshade命令-在图形中加入经济周期阴影
• Stata：图形美颜-自定义绘图模板-grstyle-palettes
• Stata绘图：多维柱状图绘制
• Stata绘图：用暂元统一改变图形中的字号
• Stata绘图全解：绘图语法-条形图-箱型图-散点图-矩阵图-直方图-点图-饼图
• Stata绘图：绘制单个变量的时序图