10  graphics 入门

不是把每个绘图函数都挨个讲一遍，也不是把它们统统归纳总结，而是比较深入地介绍一、两种图形，一、两个例子，重点阐述 Base R 的绘图特点，使用图形时，注意图形本身的作用，最终，希望读者能够达到举一反三的效果。

基础绘图系统。相比于 ggplot2 和 lattice，graphics 制作示意图是优势。

10.1 绘图基础

利用点、线等基础元素从零开始绘图。

10.1.1 plot()

本节将主要基于鸢尾花数据集介绍 R 语言基础绘图系统，该数据集最早来自埃德加·安德森，后来，被罗纳德·费希尔在介绍判别分析的论文中用到，从而，流行于机器学习社区。鸢尾花是非常漂亮的一种花，在统计和机器学习社区家喻户晓，更别提在植物界的名声。其实，远不止于此，在绘画艺术界也是如雷贯耳，印象派大师文森特·梵高画了一系列鸢尾花作品。万紫千红，但能入画的不多，故而，鸢尾花更显高雅。在生命最后的一段日子里，梵高受精神病折磨，在法国普罗旺斯的圣·雷米医院里，唯有盛开的鸢尾花陪着他，最著名的《星月夜》就是在这时候创作出来的。下面先让我们一睹鸢尾花芳容，图片来自维基百科鸢尾花词条。

(a) versicolor 杂色鸢尾

(b) setosa 山鸢尾

(c) virginica 弗吉尼亚鸢尾

图 10.1: 三种鸢尾花

鸢尾花数据集已经打包在 R 软件中，而且默认已经加载到命名空间，下面用函数 summary() 查看其概况。

summary(iris)

#>   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
#>  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
#>  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
#>  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
#>  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
#>  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
#>  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
#>        Species  
#>  setosa    :50  
#>  versicolor:50  
#>  virginica :50  
#>                 
#>                 
#> 

函数 plot() 采用公式语法可以快速作图。

plot(Sepal.Length ~ Sepal.Width, data = iris)
plot(iris$Sepal.Width, iris$Sepal.Length, panel.first = grid())

(a) 公式语法绘制散点图

(b) 带背景参考线的散点图

图 10.2: 快速作图函数 plot()

函数 plot() 是一个泛型函数，传递不同类型的参数值会调用不同的绘图方法，而不同的绘图方法的参数是不同的。当采用公式语法绘图时，会自动调用函数 plot.formula() ，此时，参数 panel.first 就不起作用。当不使用公式语法时，会调用函数 plot.default() ，此时，参数 panel.first 就起作用，利用该参数可以添加背景参考线。

10.1.2 标签

函数 plot() 的参数 xlab 、ylab 和 main 可以分别设置坐标轴横、纵标签和图主标题。

plot(
  Sepal.Length ~ Sepal.Width, data = iris, 
  xlab = "Sepal Width", ylab = "Sepal Length",
  main = "Edgar Anderson's Iris Data"
)

图 10.3: 标签

10.1.3 字体

作图函数 plot() 和 title() 都有参数 family ，设置该参数可以调整图形中的字体。下 图 10.4 的横纵坐标轴标签和图标题设为宋体，坐标轴刻度标签设为无衬线字体。

plot(Sepal.Length ~ Sepal.Width, data = iris, ann = FALSE, family = "sans")
title(
  xlab = "萼片宽度", ylab = "萼片长度",
  main = "埃德加·安德森的鸢尾花数据", family = "Noto Serif CJK SC"
)

图 10.4: 字体

10.1.4 分组

分组有两种方式，其一按照数据中的分类变量分组，其二按照一定的规则分组。而图形表达的方式可以借助颜色或图形元素的样式。

函数 plot() 的参数 col 和 pch 都可以用来分组，前者通过颜色，后者通过点的类型。简单起见，将数据集 iris 中的 Species 列传递给参数 col ，实现不同种类的鸢尾花配以不同的颜色。

plot(Sepal.Length ~ Sepal.Width, data = iris, col = Species, pch = 16)

图 10.5: 分组

下面采用一个简单规则将数据分成两组，将鸢尾花中 setosa 山毛榉类型且 Sepal.Length 萼片长度大于 5 厘米的分成一组，以红色填充散点代表这部分数据，与余下的散点形成对比，达到区分的目的。

plot(Sepal.Length ~ Sepal.Width, data = iris)
points(Sepal.Length ~ Sepal.Width, data = iris, 
  col = "#EA4335", pch = 16,
  subset = Species == "setosa" & Sepal.Length > 5
)

图 10.6: 分组

10.1.5 配色

经过探查，知道数据集 iris 中的 Species 列有三种取值。调用函数 palette() 设置一个超过 3 种颜色的调色板可以实现自定义配色。首先来看看当前调色板的颜色。

palette()

#> [1] "black"   "#DF536B" "#61D04F" "#2297E6" "#28E2E5" "#CD0BBC" "#F5C710"
#> [8] "gray62"

一共是 8 种颜色，效果预览见 图 10.7 。

图 10.7: 默认调色板

设置新的调色板也是用函数 palette() ，参数 value 设置新的颜色值向量，下面依次是红、蓝、绿、黄四种颜色。

palette(value = c("#EA4335", "#4285f4", "#34A853", "#FBBC05"))

函数 plot() 的调色板默认来自函数 palette() ，经过上面的调整，同一行绘图代码出来不同的效果，即 图 10.5 变成 图 10.8 。

plot(Sepal.Length ~ Sepal.Width, data = iris, col = Species, pch = 16)

图 10.8: 配色

10.1.6 注释

函数 text() 可以在图上任意位置添加文本或公式。下图在位置 (4,6.5) 处添加红色的文字 flower。

plot(Sepal.Length ~ Sepal.Width, data = iris)
text(x = 4, y = 6.5, labels = "flower", col = "#EA4335")

图 10.9: 注释

10.1.7 图例

函数 plot() 不会自动添加图例，需要使用函数 legend() 添加图例。

plot(Sepal.Length ~ Sepal.Width, data = iris, col = Species, pch = 16)
legend(x = "topright", title = "Species",
  legend = unique(iris$Species), box.col = NA, bg = NA,
  pch = 16, col = c("#EA4335", "#4285f4", "#34A853")
)

图 10.10: 图例

图例放置在绘图区域以外，比如右边空区域。此时，通过点和文本构造图例。

op <- par(mar = c(4, 4, 3, 6))
plot(
  Sepal.Length ~ Sepal.Width, data = iris, 
  col = Species, pch = 16, main = "Edgar Anderson's Iris Data"
)
text(x = 4.7, y = 6.75, labels = "Species", pos = 4, offset = .5, xpd = T)
points(x = 4.7, y = 6.5, pch = 16, cex = 1, col = "#EA4335", xpd = T)
text(x = 4.7, y = 6.5, labels = "setosa", pos = 4, col = "#EA4335", xpd = T)
points(x = 4.7, y = 6.3, pch = 16, cex = 1, col = "#4285f4", xpd = T)
text(x = 4.7, y = 6.3, labels = "versicolor", pos = 4, col = "#4285f4", xpd = T)
points(x = 4.7, y = 6.1, pch = 16, cex = 1, col = "#34A853", xpd = T)
text(x = 4.7, y = 6.1, labels = "virginica", pos = 4, col = "#34A853", xpd = T)
on.exit(par(op), add = TRUE)

图 10.11: 图例

在函数 plot() 内设置较宽的坐标轴范围，获得一个较宽的绘图区域，再用函数 points() 添加数据点，最后，使用函数 legend() 添加图例。

plot(
  x = c(2, 6), y = range(iris$Sepal.Length), type = "n",
  xlab = "Sepal Width", ylab = "Sepal Length",
  main = "Edgar Anderson's Iris Data"
)
points(Sepal.Length ~ Sepal.Width,
  col = Species, pch = 16, data = iris
)
legend(x = "right", title = "Species",
  legend = unique(iris$Species), box.col = NA, bg = NA,
  pch = 16, col = c("#EA4335", "#4285f4", "#34A853")
)

图 10.12: 图例

10.1.8 统计

添加分组线性回归线。按鸢尾花种类分组，线性回归模型拟合数据，抽取回归系数。首先，使用函数 split() 将数据集 iris 按变量 Species 分组拆分，得到一个列表，每个元素都是数据框。接着，调用函数 lapply() 将函数 lm() 作用到列表的每个元素上，得到一个列表，每个元素都是线性拟合对象。最后，再调函数 lapply() 将函数 coef() 应用到列表的每个元素上，得到回归模型的系数向量。

lapply(
  lapply(
    split(iris, ~Species), lm,
    formula = Sepal.Length ~ Sepal.Width
  ),
  coef
)

#> $setosa
#> (Intercept) Sepal.Width 
#>   2.6390012   0.6904897 
#> 
#> $versicolor
#> (Intercept) Sepal.Width 
#>   3.5397347   0.8650777 
#> 
#> $virginica
#> (Intercept) Sepal.Width 
#>   3.9068365   0.9015345

走到绘图这一步，往往是画什么内容比较清楚，分类数量、调色板都确定下来了。大致来说分 6 步：第一步，实现分组线性回归拟合；第二步，绘制分组散点图；第三步，添加分组回归线；第四步，添加图例并调整图例的位置；第五步，设置图形边界等绘图参数；第六步，添加背景网格线。输入线性拟合对象给函数 abline() 可以直接绘制回归线，不需要从拟合对象中提取回归系数。调用函数 par() 设置图形边界，特别是增加图形右侧边界以容纳图例，再调用函数 legend() 要设置 xpd = TRUE 以允许图例超出绘图区域。

# 分组线性拟合
iris_lm <- lapply(
  split(iris, ~Species), lm, formula = Sepal.Length ~ Sepal.Width
)
# 将分组变量和颜色映射
cols <- c("setosa" = "#EA4335",  "versicolor" = "#4285f4", "virginica" = "#34A853")
# 设置图形边界以容纳标签和图例
op <- par(mar = c(4, 4, 3, 8))
# 绘制分组散点图
plot(
  Sepal.Length ~ Sepal.Width,
  data = iris, col = Species, pch = 16,
  xlab = "Sepal Width", ylab = "Sepal Length",
  main = "Edgar Anderson's Iris Data"
)
# 添加背景参考线
grid()
# 添加回归线
for (species in c("setosa", "versicolor", "virginica")) {
  abline(iris_lm[[species]], col = cols[species], lwd = 2)
}
# 添加图例
legend(
  x = "right", title = "Species", inset = -0.4, xpd = TRUE,
  legend = unique(iris$Species), box.col = NA, bg = NA, lty = 1, lwd = 2,
  pch = 16, col = c("#EA4335", "#4285f4", "#34A853")
)
# 恢复图形参数设置
on.exit(par(op), add = TRUE)

图 10.13: 分组线性回归

10.2 绘图进阶

10.2.1 组合图形

点、线、多边形组合

x <- seq(-10, 10, length = 400)
y1 <- dnorm(x)
y2 <- dnorm(x, m = 3)
op <- par(mar = c(5, 4, 2, 1))
plot(x, y2,
  xlim = c(-3, 8), type = "n",
  xlab = quote(Z == frac(mu[1] - mu[2], sigma / sqrt(n))),
  ylab = "Density"
)
polygon(c(1.96, 1.96, x[240:400], 10),
  c(0, dnorm(1.96, m = 3), y2[240:400], 0),
  col = "grey80", lty = 0
)
lines(x, y2)
lines(x, y1)
polygon(c(-1.96, -1.96, x[161:1], -10),
  c(0, dnorm(-1.96, m = 0), y1[161:1], 0),
  col = "grey30", lty = 0
)
polygon(c(1.96, 1.96, x[240:400], 10),
  c(0, dnorm(1.96, m = 0), y1[240:400], 0),
  col = "grey30"
)
legend(x = 4.2, y = .4,
  fill = c("grey80", "grey30"),
  legend = expression(
    P(abs(Z) > 1.96, H[1]) == 0.85,
    P(abs(Z) > 1.96, H[0]) == 0.05
  ), bty = "n"
)
text(0, .2, quote(H[0]:~ ~ mu[1] == mu[2]))
text(3, .2, quote(H[1]:~ ~ mu[1] == mu[2] + delta))
on.exit(par(op), add = TRUE)

图 10.14: 正态总体下两样本均值之差的检验

10.2.2 多图布局

布局函数 layout() 和图形参数设置函数 par()

data(anscombe)
form <- sprintf("y%d ~ x%d", 1:4, 1:4)
fit <- lapply(form, lm, data = anscombe)
op <- par(mfrow = c(2, 2), mgp = c(2, 0.7, 0), 
          mar = c(3, 3, 1, 1) + 0.1, oma = c(0, 0, 2, 0))
for (i in 1:4) {
  plot(as.formula(form[i]),
    data = anscombe, col = "black",
    pch = 20, xlim = c(3, 19), ylim = c(3, 13),
    xlab = as.expression(substitute(x[i], list(i = i))),
    ylab = as.expression(substitute(y[i], list(i = i))),
    family = "sans"
  )
  abline(fit[[i]], col = "black")
  text(
    x = 7, y = 12, family = "sans",
    labels = bquote(R^2 == .(round(summary(fit[[i]])$r.squared, 3)))
  )
}
mtext("数据集的四重奏", outer = TRUE)
on.exit(par(op), add = TRUE)

图 10.15: 数据可视化很重要

10.3 图形选择

以不同的二维或三维图形可视化同一份多元数据。颜色图、透视图、等值线图和填充等值线图存在某种相似性，又有区别。

10.3.1 颜色图

\[ f(x,y) = \begin{cases} \frac{\sin(\sqrt{x^2 + y^2})}{\sqrt{x^2 + y^2}}, & (x,y) \neq (0,0)\\ 1, & (x,y) = (0,0) \end{cases} \]

y <- x <- seq(from = -8, to = 8, length.out = 51)
z <- outer(x, y, FUN = function(x, y) sin(sqrt(x^2 + y^2)) / sqrt(x^2 + y^2))
z[26, 26] <- 1

将绘图区域划分成网格，每个小网格对应一个颜色值。函数 image() 绘制颜色图

image(x = x, y = y, z = z, xlab = "$x$", ylab = "$y$")

图 10.16: 颜色图

10.3.2 透视图

函数 persp() 绘制透视图

op <- par(mar = c(0, 1, 2, 1))
persp(
  x = x, y = y, z = z, main = "二维函数的透视图",
  theta = 30, phi = 30, expand = 0.5, col = "lightblue",
  xlab = "$x$", ylab = "$y$", zlab = "$f(x,y)$"
)
on.exit(par(op), add = TRUE)

图 10.17: 透视图

10.3.3 等值线图

地理上，常用等高线图描述地形，等高线图和等值线图其实是一个意思。函数 contour() 绘制等值线图。

contour(x = x, y = y, z = z, xlab = "$x$", ylab = "$y$")

图 10.18: 等值线图

10.3.4 填充等值线图

函数 filled.contour() 绘制填充等值线图。

filled.contour(
  x = x, y = y, z = z, asp = 1,
  color.palette = hcl.colors,
  plot.title = {
    title(
      main = "二维函数的填充等值线图",
      xlab = "$x$", ylab = "$y$"
    )
  },
  plot.axes = {
    grid(col = "gray")
    axis(1, at = 2 * -4:4, labels = 2 * -4:4)
    axis(2, at = 2 * -4:4, labels = 2 * -4:4)
    points(0, 0, col = "blue", pch = 16)
  },
  key.axes = {
    axis(4, seq(-0.2, 1, length.out = 9))
  }
)

图 10.19: 填充等值线图

10.4 总结

10.4.1 plot2 包

plot2 包扩展 Base R 函数 plot() 的功能，在公式语法方面和 lattice 包很接近。另一个值得一提的 R 包是 basetheme ，用来设置 Base R 绘图主题。

library(plot2)
plot2(Sepal.Length ~ Sepal.Width | Species, data = iris, 
      palette = "Tableau 10", pch = 16)

图 10.20: plot2 包绘制分组散点图

或者使用参数 by 指定分组变量，效果和上图一样。

with(iris, {
  plot2(y = Sepal.Length, x = Sepal.Width, by = Species,
      palette = "Tableau 10", pch = 16)
})

还可以使用参数 legend 调整图例的位置，比如放置在绘图区域下方。

op <- par(mar = c(5, 4, .5, .5))
plot2(Sepal.Length ~ Sepal.Width | Species,
  data = iris, palette = "Tableau 10", pch = 16,
  legend = legend("bottom!", title = "Species of iris", bty = "o")
)
on.exit(par(op), add = TRUE)

图 10.21: plot2 包调整图例位置

还可以绘制其它类型的图形，如分组密度曲线图等。

with(iris, plot2(
  density(Sepal.Length), by = Species,
  bg = "by",   # 分组填充
  grid = TRUE, # 背景网格
  palette = "Tableau 10",
  legend = list("topright", bty = "o") # 右上角图例
))

图 10.22: plot2 包绘制密度曲线图

10.4.2 plot3D 包

虽然不提倡大量使用三维图形，但如何绘制三维图形却是生生不息的命题，以下仅是 R 语言社区的冰山一角。

• plotrix (Lemon 2006) 一个坐落于 R 的红灯区的 R 包。基于 Base R 各类绘图函数。
• scatterplot3d (Ligges 和 Mächler 2003) 基于 Base R 绘制三维散点图。
• misc3d (Feng 和 Tierney 2008) 绘制三维图形的杂项，支持通过 Base R、 tcltk 包和 rgl 包渲染图形。
• plot3D (Soetaert 2021) 依赖 misc3d 包，加强 Base R 在制作三维图形方面的能力。

举个比较新颖的一个例子，plot3D 包的函数 image2D() 绘制二维颜色图，细看又和 image() 函数不同，渲染出来的图形有三维的立体感。归根结底，很多时候束缚住自己的不是工具，而是视野和思维。以奥克兰 Maunga Whau 火山地形数据 volcano 为例。

library(plot3D)
image2D(volcano,
  shade = 0.2, rasterImage = TRUE, asp = 0.7,
  xlab = "南北方向", ylab = "东西方向",
  main = "奥克兰 Maunga Whau 地形图", clab = "高度",
  contour = FALSE, col = hcl.colors(100),
  colkey = list(
    at = 90 + 20 * 0:5, labels = 90 + 20 * 0:5,
    length = 1, width = 1
  )
)
op <- par(mar = c(1, 1.5, 0, 0))
persp3D(
  x = 1:87, y = 1:61, z = volcano, col = hcl.colors(100),
  ticktype = "detailed", colkey = FALSE, expand = 1, 
  phi = 35, theta = 125, bty = "b2", shade = TRUE,
  ltheta = 100, lphi = 45,
  xlab = "\n南北方向", ylab = "\n东西方向", zlab = "\n高度"
)
on.exit(par(op), add = TRUE)

(a) 函数 image2D() 二维颜色图

(b) 函数 persp3D() 三维透视图

图 10.23: 奥克兰火山地形图

值得一提，Python 社区的绘图模块 matplotlib 同样具有强大的绘图能力，三维图形也不在话下。不过，不同的绘图系统所采用的透视法不同，如下图所示。

from matplotlib import cm
from matplotlib.colors import LightSource
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# 设置 PGF 后端渲染图形
import matplotlib as mpl
mpl.use("pgf")
# XeLaTeX 编译图形
plt.rcParams.update({
    "text.usetex": True,
    "pgf.texsystem": "xelatex",
    "pgf.rcfonts": False,    # don't setup fonts from rc parameters
    "pgf.preamble": "\n".join([
            r"\usepackage[fontset=fandol,UTF8]{ctex}",
        ]),
})
# 读取数据
volcano = pd.read_csv("data/volcano.csv", header=None)
# DataFrame 转 Array
z = volcano.to_numpy()
# 数据行、列数
nrows, ncols = z.shape
# 线性序列
x = np.linspace(1, 87, ncols)
y = np.linspace(1, 61, nrows)
# 类似 R 语言函数 expand.grid()
xv, yv = np.meshgrid(x, y)
# 设置主题
fig, ax = plt.subplots(subplot_kw=dict(projection="3d"))
# 观察视角
ax.view_init(azim=30, elev=30)
# 设置坐标轴标签
ax.set_xlabel(r"南北方向", rotation=45)
ax.set_ylabel(r"东西方向", rotation=-15)
ax.set_zlabel(r"高度", rotation=90)
# 去掉多余的边空
fig.set_tight_layout(True)
# 光源照射的角度
ls = LightSource(270, 45)
# 自定义调色板
rgb = ls.shade(z, cmap=cm.viridis, vert_exag=0.1, blend_mode="soft")
# 三维透视图
surf = ax.plot_surface(
    xv, yv, z, rstride=1, cstride=1, facecolors=rgb,
    linewidth=0, antialiased=False, shade=False
)
# 渲染
plt.show()

图 10.24: matplotlib 绘制三维透视图

Feng, Dai, 和 Luke Tierney. 2008. 《Computing and Displaying Isosurfaces in R》. Journal of Statistical Software 28 (1). https://doi.org/10.18637/jss.v028.i01.

Lemon, Jim. 2006. 《plotrix: a package in the red light district of R》. R-News 6 (4): 8–12.

Ligges, Uwe, 和 Martin Mächler. 2003. 《scatterplot3d: An R Package for Visualizing Multivariate Data》. Journal of Statistical Software 8 (11): 1–20. https://doi.org/10.18637/jss.v008.i11.

Soetaert, Karline. 2021. plot3D: Plotting Multi-Dimensional Data. https://CRAN.R-project.org/package=plot3D.