如何使用 Stata 绘制不等宽柱状图？

如何使用 Stata 绘制不等宽柱状图？这是 Statalist 上的一个问题。How can I get a histogram with varying bin widths?。一种方法是使用 Nicholas J. Cox 编写的 eqprhistogram 命令，不过使用 twoway bar 也可以直接绘制出来。

链接：https://www.stata.com/support/faqs/graphics/histograms-with-varying-bin-widths/

最近也有培训班的小伙伴问到了这个问题，正好一块看下这个问题的解决方案。

首先我们看一下 eqprhistogram 命令的安装和使用：

安装 eqprhistogram

从 ssc 上安装即可：

ssc install eqprhistogram

使用方法很简单，例如：

use womenwage, clear
* bin(10) 表示指定每个条形代表分布中总概率的 1/10
eqprhistogram wage, bin(10) ///
    plot(kdensity wage, biweight w(5)) ///
    yla(, format(%6.2f)) yti(Density) ///
    leg(pos(2) ring(0)) ///
    xti("Wages in 1000s of dollars")

所以这里的每个柱体代表的数据是总体的十分之一。

而实际上如果你了解 twoway bar 的一些特殊功能，你可以直接绘制不等宽柱形图，例如对于 Altman (1991, 25)的数据集：

Age	Frequency
0-4	28
5-9	46
10-15	58
16	20
17	31
18-19	64
20-24	149
25-59	316
60+	103

我们把这个数据集输入 Stata：

clear all
input Age Frequency
 0  28
 5  46
10  58
16  20
17  31
18  64
20 149
25 316
60 103
80   .
end
sum Freq
* 计算密度
gen Density = Freq / (`r(sum)' * (Age[_n+1] - Age))
twoway bar Density Age, bartype(spanning) bstyle(histogram)

spanning将柱体向右拓展，这也就是为什么要指定上下界。 另外数据应该是排好序的，bstyle(histogram)选项不是必须的。如果柱体不是从 0 开始，你可能需要添加选项yscale(range(0))。

因此实际上我们也可以使用 twoway bar 绘制图 1：

use womenwage, clear
sum wage, meanonly

* 计算分位数
gen quantile = r(min) in 1
replace quantile = r(max) in 11
_pctile wage, nq(10)
ret list
forval i = 1/9 {
 replace quantile = r(r`i') in `=`i' + 1'
}

* 计算密度
gen density = 1 / (10 * (quantile[_n+1] - quantile))
label var density "Density"
sort quantile density
twoway bar density quantile, ///
 bartype(spanning) bstyle(histogram) ///
 yscale(range(0)) || ///
 kdensity wage, biweight w(5) ///
 yla(, format(%6.2f)) yti(Density) ///
    leg(pos(2) ring(0)) ///
    xti("Wages in 1000s of dollars")

那么如何绘制这样的一幅不等宽柱形图呢？

这幅图来源于：https://www.highcharts.com.cn/demo/highcharts/variwide

下面我们使用 twoway bar 来绘制这幅图。

首先输入数据：

clear
input str15 country cost gdp
"挪威" 50.2 335504
"丹麦" 42 277339
"比利时" 39.2 421611
"瑞典" 38 462057
"法国" 35.6 2228857
"荷兰" 34.3 702641
"芬兰" 33.2 215615
"德国" 33.0 3144050
"奥地利" 32.7 349344
"爱尔兰" 30.4 275567
"意大利" 27.8 1672438
"英国" 26.7 2366911
"西班牙" 21.3 1113851
"希腊" 14.2 175887
"葡萄牙" 13.7 184933
"捷克共和国" 10.2 176564
"波兰" 8.6 424269
"罗马尼亚" 5.5 169578
end

根据 GDP 分组生成 x 轴变量：

* 生成累加值
gen cumsum = sum(gdp)
set obs 19
gen cumsum2 = .

forval i = 2/`=_N' {
    replace cumsum2 = cumsum[`i' - 1] in `i'
}
replace cumsum2 = 0 if missing(cumsum2)

然后就可以绘制一幅简单的不等宽柱形图了：

twoway bar cost cumsum2, bartype(spanning) ///
    bstyle(histogram) yla(, format(%6.2f))

每一个柱形图的中心值为：

gen x = .
forval i = 2/`=_N' {
    replace x = (cumsum2[_n + 1] + cumsum2[_n]) / 2
}

再生成一个 color 变量：

配色网站：https://tidyfriday.cn/colors/

input str11 color
"230 75 53"
"77 187 213"
"0 160 135"
"60 84 136"
"243 155 127"
"132 145 180"
"145 209 194"
"220 0 0"
"126 97 72"
"176 156 133"
"0 160 135"
"60 84 136"
"243 155 127"
"132 145 180"
"145 209 194"
"220 0 0"
"126 97 72"
"176 156 133"
end

散点标签：

tostring cost, gen(costlab) format(%6.0f) force
replace costlab = costlab + " €/h"

循环生成绘图语句：

local cmd = "tw"
local leg = ""
local xtext = ""
forval i = 1/18 {
    * 绘图语句
    local cmd = `"`cmd' (bar cost cumsum2 in `i'/`=`i'+1', bartype(spanning) bstyle(histogram) fc("`=color[`i']'") lc("`=color[`i']'"))"'
    * 图例
    local leg = `"`leg' `i' "`=country[`i']'" "'
    * x 轴标签
    if "`=country[`i']'" != "葡萄牙"{
        local xtext = `" `xtext' `=x[`i']' "`=country[`i']'" "'
    }
}

`cmd', xla(`xtext', ang(90) labsize(*0.9) nogrid) ///
    xti("") yti("人工成本(€/h)") || ///
    sc cost x if gdp >= 1113851, ///
    m(i) mla(costlab) mlabpos(12) mlabcolor(black) ///
    yla(0(8)56) ///
    ti("2016 年欧洲各国人工成本") ///
    subti("绘制：微信公众号 RStata｜柱体的宽度表示各国 GDP 的大小") ///
    caption("数据来源：eurostat" "http://ec.europa.eu/eurostat/web/") ///
    leg(off)

这里我虽然生成了图例标签，但是并没有使用。

点击这里跳转到 RStata 短书平台获取附件：如何使用 Stata 绘制不等宽柱状图？