R语言之基础绘图

文章和代码已经归档至【Github仓库：https://github.com/timerring/dive-into-AI 】或者公众号【AIShareLab】回复 R语言 也可获取。

R 的基础绘图系统由 Ross Ihaka 编写，功能非常强大，主要由 graphics 包和 grDevices 包组成，它们在启动 R 时会自动加载。基础绘图系统中有两类函数，一类是高水平作图函数，另一类是低水平作图函数。

所谓高水平作图函数是用于直接产生图形的函数，包括 plot( )、hist( )、boxplot( )和 pairs( )等。低水平作图函数是用于在高水平作图函数所绘图形的基础上添加新的图形或元素的函数，包括 points( )、lines( )、text( )、title( )、legend( )和 axis( )等。

1. 函数 plot( )

函数 plot( ) 是一个泛型函数，对于不同类型的数据，它可以绘制出不同的图形。例如，对于数值型数据，它可以绘制出散点图；对于分类数据，它可以绘制出箱线图；对于一些统计模型，它可以绘制出相应的图形，比如对于生存分析，它可以绘制出生存曲线。因此，函数 plot( ) 的使用频率非常高，建议可以打开它的帮助文档查看其各种常用参数的用法。

下面创建一个示例数据，表示某病病人对 2 种药物（drugA 和 drugB）、5 个剂量（dose）水平上的响应情况。

dose <- c(20, 30, 40, 45, 60)
drugA <- c(16, 20, 27, 40, 60)
drugB <- c(15, 18, 25, 31, 40)

用上面的数据绘制药物 A 的剂量和响应关系的图形：

plot(dose, drugA)
plot(dose, drugA, type = "b")

上面的命令创建了两幅图，函数 plot( ) 里的参数 type 默认为“p”（代表点），所以得到的第一张图是散点图。在第二行命令里，参数 type 改为了“b”（代表点和线），所以得到的第二张图是点线图。

函数 plot( ) 用于新建一幅图形，我们还可以用低水平作图函数，例如 lines( )、legend( )等，在一幅现有图形上添加新的图形元素。例如：

# 为了比较两种药物不同剂量下的响应情况，我们在一幅图上展示两个点线图，并用不同类型的线（lty）和不同特征的点（pch）加以区分。
plot(dose, drugA, type = "b", lty = 1, pch = 15)
lines(dose, drugB, type = "b", lty = 2, pch = 17)
# 为了增强可读性，还添加了图例（legend）。
# 需要注意的是，函数 legend( )里面点和线的属性必须与前面函数 plot( )和 lines( )中设置的属性一致。
legend("topleft", title = "Drug Type",
       legend = c("A", "B"), 
       lty = c(1, 2), 
       pch = c(15, 17))

2.直方图和密度曲线图

直方图（histogram）是用于展示连续型变量分布的最常用的工具，它本质上是对密度函数的一种估计。直方图和密度曲线图一般用于探索分布，很少用于报告结果。函数 hist( )可用于绘制直方图。

数据集 anorexia 位于 MASS 包中，来自一项关于年轻女性厌食症患者体重变化的研究。该数据集包含 72 例观察对象、3 个变量，其中变量 Treat（治疗方式）是一个包含 3 个水平的因子，变量 Prewt 和 Postwt 均为数值型，分别表示治疗前后的体重（单位：lb）。下面绘制变量 Prewt 的直方图，代码如下：

library(MASS)
data(anorexia)
str(anorexia)

attach(anorexia)
hist(Prewt)

上图给出了变量 Prewt 的频数分布，由于函数 hist( )中没有设置任何参数，图中使用了默认的组距、坐标轴标签和标题等。需要注意的是，直方图的形状受到组距的影响，有时我们需要尝试设定参数 breaks 的不同的值以得到合适的图形。函数 hist( )的输出结果中包含一些计算返回值，这些值可用于进一步地作图或者分析，例如为区间划分端点、频数（或密度）、区间中点等。

密度曲线为数据的分布提供了一种更为平滑的描述，绘制密度曲线的方法为：

plot(density(Prewt))

从上图可以看出，变量 Prewt 的分布是单峰的，基本是对称的。我们还可以在一幅直方图上添加一条密度曲线和轴须图。此时，需要在函数 hist( )里面设定参数 freq 为 FALSE，即把纵坐标换成频率，否则将会几乎看不到密度曲线。参数 las（或者 labels） 设为 1 是为了将纵轴的刻度标签横向显示。

library("showtext") # R 数据分析镜像的中文支持不太好，需要借助 showtext 包
showtext_auto() # 自动支持中文
# 使用红色填充了条形，添加了信息量更大的坐标轴标签和标题，还通过设置参数 las 为 1 把纵轴的刻度标签换成了横向显示。
hist(Prewt, freq = FALSE, col = "red",
     xlab = "体重(lbs)", 
     main = "治疗前体重分布直方图",
     las = 1)
# 然后使用函数 lines( )在直方图上叠加了一条蓝色的、两倍于默认线条宽度的密度曲线。
lines(density(Prewt), col = "blue", lwd = 2)
# 最后使用函数 rug( )在横轴上添加了轴须图，以展示数据分布的密集趋势。
rug(Prewt)
detach(anorexia)

3.条形图

条形图（bar chart）在医学科技论文中经常用到，它通过垂直的或水平的矩形展示分类变量的频数分布。函数 barplot( ) 可用于绘制条形图。

下面以 vcd 包里的 Arthritis 数据集为例介绍函数 barplot( )的用法。该数据集来自一项关于治疗类风湿性关节炎新方法的成组对照双盲临床试验研究。其中的反应变量 Improved 记录了每位接受药物治疗（Treated，41 例）或安慰剂（Placebo，43 例）的患者的治疗效果，分为 3 个级别（None、Some、Marked）。

library(vcd)
data(Arthritis)
attach(Arthritis)
counts <- table(Improved)
counts
# Improved
#   None   Some Marked 
#     42     14     28

函数 table( )用于生成分类变量的频数统计表。从上面的输出可以看到，有 28 位患者有了明显改善、14 人有部分改善，而有 42 人没有改善。条形图可以用于展示这个频数分布，如下图所示：

barplot(counts, xlab = "Improvement", ylab = "Freqency", las = 1)

函数 barplot( )还可以用于展示二维列联表的数据。下图绘制了一幅分组条形图，并添加了颜色和图例，代码如下：

counts <- table(Improved, Treatment)
barplot(counts, 
        col = c("red", "yellow", "green"),
        xlab = "Improvement", 
        ylab = "Freqency", 
        beside = TRUE, las = 1)
legend("top", legend = rownames(counts), 
       fill = c("red", "yellow", "green"))

条形图有时候也可以用于展示不同分类下的均值、中位数、标准差、置信区间等。用基本包里的函数可以实现这个功能，但是需要很多个步骤。而 epiDisplay 包里的函数 aggregate.plot( )可以简化这个过程。

下面的代码以数据集 anorexia 为例绘制了不同治疗方式下治疗后体重的均值条形图，结果如下图所示。

library(epiDisplay)
aggregate.plot(anorexia$Postwt, by = list(anorexia$Treat), 
               error = "sd", legend = FALSE, 
               bar.col = c("red", "yellow", "green"),
               ylim = c(0,100), las = 1,
               main = "")

上面的误差棒表示的是标准差，我们可以通过改变函数 aggregate.plot( )里面的参数 error 设置显示标准误或置信区间。

4. 饼图

饼图（pie chart）可用于展示分类数据的占比情况。例如，下面的代码绘制的饼图展示了某医院一周内急诊入院的疾病类型分布。

percent <- c(5.8, 27.0, 0.5, 20.8, 12.8, 33.1)
disease <- c("上感", "中风", "外伤", "昏厥", "食物中毒", "其他")
lbs <- paste0(disease, percent, "%")
pie(percent, labels = lbs, col = rainbow(6))

多数统计学家不建议使用饼图，他们更推荐用条形图或点图代替饼图，因为人们对长度的判断比对面积的判断更精确。因此，基本包的函数 pie( )绘制饼图的选项有限。

不过，一些捐赠包扩展了 R 绘制饼图的功能，例如 plotrix 包。该包提供的函数 pie3D( )可以绘制三维饼图，另一个函数 fan.plot( )可以绘制功能与饼图相似的扇形图，感兴趣的读者可以安装该包并查看其帮助文档。

5. 箱线图和小提琴图

箱线图（box plot）又称箱须图（box-whisker plot），常用于展示数据的大致分布特征，也用于探索异常值和离群点。函数 boxplot( )可用于绘制箱线图。

下面用箱线图展示数据集 anorexia 里体重前后变化的分布。

anorexia$wt.change <- anorexia$Postwt - anorexia$Prewt
boxplot(anorexia$wt.change, ylab = "Weight change (lbs)", las = 1)

为了让读者更好地理解箱线图各部分的含义，在下图中添加了手工标注。如果数据是对称分布，中位数（Median）应该位于上四分位数（Upper quartile）和下四分位数（Lower quartile）的中间，即箱线图的方盒关于中位线对称。在上边缘（Upper hinge）和下边缘（Lower hinge）以外的值通常被认为是异常值。

fivenum(anorexia$wt.change)

anorexia$wt.change <- anorexia$Postwt - anorexia$Prewt
b <- boxplot(anorexia$wt.change, ylab = "Weight change (lbs)", las = 1)
# text(x= 1, y=1:5, labels= c("some","more","red text"))
text(1.2, 21.5, "Upper hinge")
text(1.13, 15.5, "←—— Whisker")
text(1.31, 9.2, "Upper quantile")
text(1.26, 1.65, "Median")
text(1.31, -2.45, "Upper quantile")
text(1.13, -7, "←—— Whisker")
text(1.2, -12.2, "Upper hinge")

平行排列的箱线图可以用于比较在某个分类变量各个类别下某指标的分布。例如，要比较不同治疗方式下体重变化的情况，可以使用下面的命令：

boxplot(wt.change ~ Treat, data = anorexia,
        ylab = "Weight change (lbs)", las = 1)

函数 boxplot( )的第一个参数输入的是一个公式。R 里公式一般用符号 ~ 连接变量，~ 左边可以看作因变量， ~ 右边可以看作自变量。从下图（a）可以看出，“FT”（family treatment）组体重的改变量高于其他两组。但是，差异的显著性需要进一步的显著性检验才能确定。

小提琴图（violin plot）可以看作是箱线图和密度图的结合。vioplot 包里的函数 vioplot( )可用于绘制小提琴图，使用前请先安装并加载该包。例如，上图可以换成小提琴图展示

options(warn=-1) # 清爽显示
library(vioplot)
vioplot(wt.change ~ Treat, data = anorexia, 
        ylab = "Weight change (lbs)",
        col = "gold", las = 1)

6. 克利夫兰点图

克利夫兰点图（Cleveland dot plot）本质上也是散点图，它通过点的位置展示数据的大小，是一种在简单水平刻度上绘制大量有标签的值的方法，其功能与条形图类似，但强调数据的排序以及相互之间的差距。

函数 dotchart( )可用于绘制克利夫兰点图。datasets 包里的数据集 VADeaths 是 1940 年美国弗吉尼亚州城市和农村不同年龄段人群的死亡率（以 ‰ 表示）情况。

VADeaths
dotchart(VADeaths)
dotchart(t(VADeaths),pch = 19)

从上图可以看出，死亡率随着年龄的升高而升高；在同一年龄段，农村地区的死亡率均高于城市地区；在同一年龄段同一地区，男性的死亡率均高于女性。

7. 导出图形

如果想要把图形保存下来，可以通过图形用户界面和代码两种方式。在 RStudio 右下方的“Plots”下，单击“Export”，选择“Save as Image”或“Save as PDF”，可以把图形保存在指定的文件夹下。我们还可以选择“Copy to Clipboard”把图形直接复制到 Word 或 PowerPoint 文档。需要注意的是，这种方式保存的图形与 RStudio 图形窗口的尺寸有关，即不同大小的窗口得出的图形会有差异（在 ModelWhale 中，可以右键单击图片，直接另存为）。

如果想把图形保存下来用于报告或论文中，笔者建议使用代码的方式，将绘图语句放置在开启目标图形设备的语句和关闭目标图形设备的语句之间即可。例如，下面的代码会把图形保存到当前工作目录中并命名“mygraph.pdf”：

pdf("mygraph.pdf")
boxplot(wt.change ~ Treat, 
        data = anorexia, 
        ylab = "Weight change (lbs)",
        las = 1)
dev.off() # work 下可以看到该 pdf

除了函数 pdf( )，我们还可以使用函数 png( )、jpeg( )、tiff( )和 postscript( )等将图形保存为其他格式。

bmp、png 和 jpeg 格式的图形文件都是非矢量格式，容易受到分辨率的影响，但它们占用的空间很小，适合运用于 Word 和PowerPoint 文档中；ps 格式的图形文件是矢量格式文件，它与分辨率无关，适合运用于排版印刷；而 tiff（或 tif）格式的图形文件可以支持很多色彩系统，而且独立于操作系统，在各类出版物中运用得最为广泛。例如：

tiff(filename = "mygraph.tiff",
     width = 15, height = 12, units = "cm", res = 300)
boxplot(wt.change ~ Treat, data = anorexia, ylab = "Weight change (lbs)")
dev.off() # work 下可以看到该 tiff

上述命令生成了一个名为“mygraph.tiff”的图形文件，参数 width 和 height 分别用于设置图形的宽度和高度，参数 units 用于设置宽度和高度的单位，参数 res 用于设置分辨率，这里设为了大多数出版物要求的最低值 300。

小结

其他一些专门的图形，例如散点图矩阵、相关图、正态 QQ 图、生存曲线、聚类图、碎石图、ROC 曲线和 Meta 分析森林图等。在 R 的应用中，可视化是一个非常活跃的领域，新的包层出不穷。网站 The R Graph Gallery 收集了各种新颖的图形以及相应的示例代码，值得对可视化感兴趣的读者关注。