之前给大家介绍了使用 R 语言绘制相关系数热力图的方法，今天我们再来学习下 Stata 如何绘制相关系数热力图。这里我总结了 3 种方法，推荐掌握第一种。
为了让大家更好的理解下面的代码，欢迎各位培训班会员参加明晚 8 点的直播课：「使用 Stata 绘制相关系数热力图！」
1. 散点图法
从另外一个角度看，热力图其实就是散点图，所以我们可以采用绘制散点图的方式来绘制热力图。
首先加载我们的示例数据 mtcars.dta（R 语言的 mtcars 数据集，保存成 mtcars.dta 即可）：
计算所有变量的相关系数矩阵和 p 值：
pwcorr, sig
* 计算的结果会被存储到返回值中：
ret list
mat a = r(C)
mat b = r(sig)
|
例如 a 是这样的：
mat list a
*> symmetric a[11,11]
*> mpg cyl disp hp drat wt qsec vs am gear carb
*> mpg 1
*> cyl -.85216196 1
*> disp -.84755138 .90203287 1
*> hp -.77616837 .83244745 .79094859 1
*> drat .68117191 -.69993811 -.71021393 -.44875912 1
*> wt -.86765938 .78249579 .88797992 .65874789 -.71244065 1
*> qsec .41868403 -.59124207 -.43369788 -.70822339 .09120476 -.17471588 1
*> vs .66403892 -.8108118 -.71041589 -.72309674 .44027846 -.55491568 .74453544 1
*> am .59983243 -.52260705 -.59122704 -.24320426 .71271113 -.69249526 -.22986086 .16834512 1
*> gear .48028476 -.4926866 -.5555692 -.12570426 .69961013 -.583287 -.21268223 .20602335 .79405876 1
*> carb -.55092507 .52698829 .39497686 .74981247 -.0907898 .42760594 -.65624923 -.56960714 .05753435 .27407284 1
|
调用矩阵 a 里面的元素可以采用下面的方法：
di a["mpg", "cyl"]
*> -.85216196
|
下面我们生成一个数据集，这个数据集有三个变量，varname1、varname2 以及他们俩的相关系数 corr，由于我们这里的矩阵是 11x11 的，所以这个数据集应该是 121 个观测值，可以通过如下代码生成：
clear
input str10 varname1
"mpg"
"cyl"
"disp"
"hp"
"drat"
"wt"
"qsec"
"vs"
"am"
"gear"
"carb"
end
save mydata, replace
ren varname1 varname2
cross using mydata
gen corr = .
gen p = .
forval i = 1/`=_N' {
replace corr = a["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
replace p = b["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
}
save mydata, replace
|
但是如果变量很多的话，input 可能就比较累了，可以采用下面的方法：
use mtcars.dta, clear
tempfile tmp
file open myfile using `tmp', write
file write myfile "varname1" _n
foreach i of varlist _all {
file write myfile "`i'" _n
}
file close myfile
import delimited using `tmp', varnames(1) clear
save mydata, replace
ren varname1 varname2
cross using mydata
gen corr = .
gen p = .
forval i = 1/`=_N' {
replace corr = a["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
replace p = b["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
}
save mydata, replace
|
然后我们就可以绘图了，不过 Stata 绘图命令不接受字符串变量，所以我们需要把 varname1 和 varname2 转成因子变量，使用 sencode 命令可以把一个字符串变量按照另外一个变量的大小顺序转换成有序因子变量：
* 安装 sencode
* ssc install sencode
sencode varname1, gsort(corr) gen(v1)
sencode varname2, gsort(corr) gen(v2)
|
然后我们再生成一个变量来标志显著性：
gen star = "" if missing(p)
replace star = "***" if p <= 0.001
replace star = "**" if inrange(p, 0.001, 0.01)
replace star = "*" if inrange(p, 0.01, 0.1)
tostring corr, format(%6.2f) gen(corr2) force
|
热力图通常用色块的颜色来表示除了 x y 的第三个变量，这里就是相关系数了，但是 Stata 不支持直接把散点的颜色用某个变量进行映射，所以我们需要分多个图层绘制，例如分 9 个，我们先生成分层变量：
egen group = cut(corr), group(9) label
|
然后就可以分 9 层绘图了：
配色的选择可以参考这个网站：https://tidyfriday.cn/colors
tw ///
sc v1 v2 if group == 0, m(S) msize(*5.2) mc("178 24 43") || ///
sc v1 v2 if group == 1, m(S) msize(*5.2) mc("214 96 77") || ///
sc v1 v2 if group == 2, m(S) msize(*5.2) mc("244 165 130") || ///
sc v1 v2 if group == 3, m(S) msize(*5.2) mc("253 219 199") || ///
sc v1 v2 if group == 4, m(S) msize(*5.2) mc("247 247 247") || ///
sc v1 v2 if group == 5, m(S) msize(*5.2) mc("209 229 240") || ///
sc v1 v2 if group == 6, m(S) msize(*5.2) mc("146 197 222") || ///
sc v1 v2 if group == 7, m(S) msize(*5.2) mc("67 147 195") || ///
sc v1 v2 if group == 8, m(S) msize(*5.2) mc("33 102 172") || ///
sc v1 v2, m(i) mlab(corr2) mlabpos(12) mlabc(balck) mlabsize(*0.8) || ///
sc v1 v2, m(i) mlab(star) mlabpos(6) mlabc(balck) mlabsize(*0.8) aspectr(1) ///
leg(off) xla(1(1)11, value nogrid notick labgap(*5)) ///
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline) ///
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
|
| 分图层绘制 |
 |
这样绘图会很麻烦，如果我们需要更多图层就更麻烦了，所以这个时候我们可以使用循环生成绘图语句：
* 使用循环生成绘图语句
local cmd = "tw"
local j = 0
foreach c in "197 27 125" "222 119 174" "241 182 218" "253 224 239" "247 247 247" "230 245 208" "184 225 134" "127 188 65" "77 146 33" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(S) msize(*5.2) mc("`c'"))"'
local j = `j' + 1
}
`cmd' || ///
sc v1 v2, m(i) mlab(corr2) mlabpos(12) mlabc(balck) mlabsize(*0.8) || ///
sc v1 v2, m(i) mlab(star) mlabpos(6) mlabc(balck) mlabsize(*0.8) aspectr(1) ///
leg(off) xla(1(1)11, value nogrid notick labgap(*5)) ///
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline) ///
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
|
| 发散配色1 |
 |
那这个时候我们就可以很容易的更换颜色了！
local cmd = "tw"
local j = 0
foreach c in "140 81 10" "191 129 45" "223 194 125" "246 232 195" "245 245 245" "199 234 229" "128 205 193" "53 151 143" "1 102 94" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(S) msize(*5.2) mc("`c'"))"'
local j = `j' + 1
}
`cmd' || ///
sc v1 v2, m(i) mlab(corr2) mlabpos(12) mlabc(balck) mlabsize(*0.8) || ///
sc v1 v2, m(i) mlab(star) mlabpos(6) mlabc(balck) mlabsize(*0.8) aspectr(1) ///
leg(off) xla(1(1)11, value nogrid notick labgap(*5)) ///
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline) ///
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
|
| 发散配色2 |
 |
也可以使用渐变色：
local cmd = "tw"
local j = 0
foreach c in "224 242 241" "178 223 219" "128 203 196" "77 182 172" "38 166 154" "0 150 136" "0 137 123" "0 121 107" "0 105 92" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(S) msize(*6) mc("`c'"))"'
local j = `j' + 1
}
`cmd' || ///
sc v1 v2, m(i) mlab(corr2) mlabpos(12) mlabc(balck) mlabsize(*0.8) || ///
sc v1 v2, m(i) mlab(star) mlabpos(6) mlabc(balck) mlabsize(*0.8) aspectr(1) ///
leg(off) xla(1(1)11, value nogrid notick labgap(*5)) ///
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline) ///
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
|
| 渐变配色 |
 |
也可以很容易的更改色块的形状：
local cmd = "tw"
local j = 0
foreach c in "140 81 10" "191 129 45" "223 194 125" "246 232 195" "245 245 245" "199 234 229" "128 205 193" "53 151 143" "1 102 94" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(O) msize(*5.2) mc("`c'"))"'
local j = `j' + 1
}
`cmd' || ///
sc v1 v2, m(i) mlab(corr2) mlabpos(12) mlabc(balck) mlabsize(*0.8) || ///
sc v1 v2, m(i) mlab(star) mlabpos(6) mlabc(balck) mlabsize(*0.8) aspectr(1) ///
leg(off) xla(1(1)11, value nogrid notick labgap(*5)) ///
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline) ///
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
|
| 圆形热力图 |
 |
实际上还可以下面这样玩（Windows 上未经测试）：
从附件中安装 ChineseZodiac 字体（双击字体文件即可安装），Windows 电脑应该是在电脑字体设置里面找到 ChineseZodiac 字体的名字然后替换下面代码里面的 “ChineseZodiac”（也就是说字体的名字可能不同）。
cap drop group
egen group = cut(corr), group(12) label
cap drop tmp
ascii
gen tmp = ""
local j = 0
forval i = 0/11 {
local j = `i' + 97
replace tmp = `"{fontface "ChineseZodiac": `=char(`j')'}"' if group == `i'
}
local cmd = "tw"
local j = 0
foreach c in "84 48 5" "140 81 10" "191 129 45" "223 194 125" "246 232 195" "245 245 245" "199 234 229" "128 205 193" "53 151 143" "1 102 94" "0 60 48" "black" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(i) mlabsize(*2.5) mlab(tmp) mlabc("`c'") mlabpos(0))"'
local j = `j' + 1
}
`cmd', aspectr(1)
leg(off) xla(1(1)11, value nogrid notick labgap(*5))
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline)
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12)
ti("mtcars 数据集相关系数矩阵", size(*1.2))
subti("数据来源：datasets 包" " ", size(*0.8))
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
graphr(margin(5 5 5 5))
|
| 十二生肖 |
 |
还可以这样玩（仅限 Mac）：
local j = 0
foreach i in 🍣 🍤 🍟 🍞 🍝 🍡 🍩 🍪 🍫 🍭 🍧 🍯 {
replace tmp = "`i'" if group == `j'
local j = `j' + 1
}
local cmd = "tw"
local j = 0
foreach c in "84 48 5" "140 81 10" "191 129 45" "223 194 125" "246 232 195" "245 245 245" "199 234 229" "128 205 193" "53 151 143" "1 102 94" "0 60 48" "black" {
local cmd = `"`cmd' (sc v1 v2 if group == `j', m(i) mlabsize(*2.5) mlab(tmp) mlabc("`c'") mlabpos(0))"'
local j = `j' + 1
}
`cmd', aspectr(1)
leg(off) xla(1(1)11, value nogrid notick labgap(*5))
yla(1(1)11, value nogrid notick labgap(*5)) ysc(noline)
xsc(noline) xti("") yti("") scheme(plotplain) xsize(10) ysize(12)
ti("mtcars 数据集相关系数矩阵", size(*1.2))
subti("数据来源：datasets 包" " ", size(*0.8))
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
graphr(margin(5 5 5 5))
|
| Mac 表情包 |
 |
2. 使用 plotmatrix 绘制
plotmatrix 是个外部命令：
* 安装
ssc install plotmatrix
|
用法也很简单：
* 计算所有变量的相关系数矩阵和 p 值
pwcorr, sig
mat a = r(C)
plotmatrix, mat(a) leg(off) formatcells(%5.2f) ///
freq maxticks(11) ///
ti("mtcars 数据集相关系数矩阵", size(*1.2)) ///
subti("数据来源：datasets 包" " ", size(*0.8)) ///
caption("绘制：微信公众号 RStata", size(*0.6) pos(5)) ///
scheme(plotplain) xsize(10) ysize(12) ///
ysc(noline) xsc(noline) xti("") yti("") ///
xla(, value nogrid notick labgap(*5)) ///
yla(, value nogrid notick labgap(*5)) ysc(noline)
|
| plotmatrix 画法 |
 |
不过标签不能自定义。
3. 作为等高图绘制
这种方法绘制前我们也需要生成一个数据集（和第一种方法一样）：
use mtcars.dta, clear
pwcorr, sig
ret list
mat a = r(C)
mat b = r(sig)
clear
input str10 varname1
"mpg"
"cyl"
"disp"
"hp"
"drat"
"wt"
"qsec"
"vs"
"am"
"gear"
"carb"
end
save mydata, replace
ren varname1 varname2
cross using mydata
gen corr = .
gen p = .
forval i = 1/`=_N' {
replace corr = a["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
replace p = b["`=varname1[`i']'", "`=varname2[`i']'"] in `i'
}
save mydata, replace
use mydata, clear
sencode varname1, gsort(corr) gen(v1)
sencode varname2, gsort(corr) gen(v2)
tw contour corr v1 v2, interp(none) heatmap
ccolors("103 0 31" "178 24 43" "214 96 77" "244 165 130" "253 219 199" "247 247 247" "209 229 240" "146 197 222" "67 147 195" "33 102 172" "5 48 97")
aspectr(1) levels(11)
xla(1(1)11, value nogrid notick labgap(*5))
yla(1(1)11, value nogrid notick labgap(*5))
ysc(noline) xsc(noline)
xti("") yti("") scheme(plotplain)
xsize(12) ysize(10)
ti("mtcars 数据集相关系数矩阵", size(*1.2))
subti("数据来源：datasets 包" " ", size(*0.8))
caption("绘制：微信公众号 RStata", size(*0.6) pos(5))
cleg(width(3)) zti("相关系数", size(*0.8))
zla(-1(0.5)1, format(%6.1f) notick)
|
| tw contour 画法 |
 |
4. 总结
可以看出后面两种方法虽然看起来很简单，但是都不如第一种方法灵活（想怎么画就怎么画），推荐掌握第一种。
点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制相关系数热力图
评论