# 第 4 章 图库

— 柯南·道尔《血字的研究》

## 4.1 直方图

$$$f(x)=F'(x)=\lim_{h\rightarrow0}\frac{F(x+h)-F(x)}{h} \tag{4.1}$$$

$$$\hat{F}_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{I}(X_{i}\leq x) \tag{4.2}$$$

$$$\hat{f}_{n}(x)=\lim_{h\rightarrow0}\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{I}(x<X_{i}\leq x+h)}{h} \tag{4.3}$$$

$$$b_{1}<b_{2}<\cdots<b_{j}<b_{j+1}<\cdots;\;b_{j+1}-b_{j}=h,\,j=1,2,\cdots \tag{4.4}$$$

$$$\hat{f}_{n}(x)=\frac{1}{nh}\sum_{i=1}^{n}\mathbf{I}(b_{j}<X_{i}\leq b_{j+1});\;x\in(b_{j},b_{j+1}] \tag{4.5}$$$

usage(hist.default)
## ## Default S3 method:
## hist(x, breaks = "Sturges", freq = NULL, probability = !freq,
##   include.lowest = TRUE, right = TRUE, density = NULL, angle = 45,
##   col = "lightgray", border = NULL, main = paste("Histogram of", xname),
##   xlim = range(breaks), ylim = NULL, xlab = xname, ylab, axes = TRUE,
##   plot = TRUE, labels = FALSE, nclass = NULL, warn.unused = TRUE, ...)

par(mfrow = c(2, 2), mar = c(2, 3, 2, .5), mgp = c(2, .5, 0))
data(geyser, package = "MASS")
hist(geyser$waiting, main = "(1) freq = TRUE", xlab = "waiting") hist(geyser$waiting, freq = FALSE, xlab = "waiting", main = "(2) freq = FALSE")
hist(geyser$waiting, breaks = 5, density = 10, xlab = "waiting", main = "(3) breaks = 5") hist(geyser$waiting, breaks = 40, col = "red", xlab = "waiting", main = "(4) breaks = 40")

library(ggplot2)
library(cowplot)
p <- ggplot(aes(waiting), data = geyser)
p1 <- p + geom_histogram(breaks = seq(40, 110, by = 5))
p2 <- p + geom_histogram(breaks = seq(40, 110, by = 5), aes(y = ..density..))
p3 <- p + geom_histogram(breaks = seq(40, 110, by = 10))
p4 <- p + geom_histogram(breaks = seq(42, 108, by = 2), fill = "red", color = "black")
plot_grid(p1, p2, p3, p4, labels = c(
"(1) freq = TRUE",
"(2) freq = FALSE",
"(3) breaks = 5",
"(4) breaks = 40"
), ncol = 2)

demo("hist_geyser", package = "MSG")
df <- data.frame(
x = seq(40, 110, 5), y = 0,
xend = seq(40, 110, 5), yend = ht
)
p2 + geom_density(fill = "lightgray", color = "black") +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend),
data = df, lty = 3
)

## 4.2 茎叶图

R 中茎叶图的函数为 stem()，其用法为：

usage(stem)
## stem(x, scale = 1, width = 80, atom = 1e-08)

head(islands, 10)
##       Africa   Antarctica         Asia    Australia Axel Heiberg       Baffin
##        11506         5500        16988         2968           16          184
##        Banks       Borneo      Britain      Celebes
##           23          280           84           73

stem(islands)
##
##   The decimal point is 3 digit(s) to the right of the |
##
##    0 | 00000000000000000000000000000111111222338
##    2 | 07
##    4 | 5
##    6 | 8
##    8 | 4
##   10 | 5
##   12 |
##   14 |
##   16 | 0
stem(islands, width = 20)
##
##   The decimal point is 3 digit(s) to the right of the |
##
##    0 | 00000000+21
##    2 | 07
##    4 | 5
##    6 | 8
##    8 | 4
##   10 | 5
##   12 |
##   14 |
##   16 | 0
# 可以增大窗宽 stem(islands, scale = 2) 看看效果

# 去掉陆地名称以便显示数据
unname(sort(round(islands / 1000, 1)))
##  [1]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
## [16]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.1
## [31]  0.1  0.1  0.1  0.1  0.1  0.2  0.2  0.2  0.3  0.3  0.8  3.0  3.7  5.5  6.8
## [46]  9.4 11.5 17.0

# 均值 lambda 为 10 的泊松分布随机数
sort(x <- rpois(80, lambda = 10))
##  [1]  5  5  6  6  6  6  6  6  7  7  7  7  7  7  7  7  7  8  8  8  8  8  8  8  8
## [26]  8  8  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9 10 10 10 10 10 10 10 10
## [51] 10 11 11 11 11 11 11 12 12 12 12 12 12 12 13 13 13 13 13 13 13 14 14 14 14
## [76] 14 14 15 16 16
stem(x, scale = 2)
##
##   The decimal point is at the |
##
##    5 | 00
##    6 | 000000
##    7 | 000000000
##    8 | 0000000000
##    9 | 000000000000000
##   10 | 000000000
##   11 | 000000
##   12 | 0000000
##   13 | 0000000
##   14 | 000000
##   15 | 0
##   16 | 00

## 4.3 箱线图

R 中相应的函数为 boxplot()，其用法如下：

# 默认用法
usage(boxplot.default)
## ## Default S3 method:
## boxplot(x, ..., range = 1.5, width = NULL, varwidth = FALSE, notch = FALSE,
##   outline = TRUE, names, plot = TRUE, border = par("fg"), col = "lightgray",
##   log = "", pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5),
##   ann = !add, horizontal = FALSE, add = FALSE, at = NULL)
# 公式用法
usage(graphics:::boxplot.formula)
## ## S3 method for class 'formula'
## boxplot(formula, data = NULL, ..., subset, na.action = NULL,
##   xlab = mklab(y_var = horizontal), ylab = mklab(y_var = !horizontal),
##   add = FALSE, ann = !add, horizontal = FALSE, drop = FALSE, sep = ".",
##   lex.order = FALSE)

boxplot(count ~ spray, data = InsectSprays,
col = "lightgray", horizontal = TRUE, pch = 4)
ggplot(aes(y = count, x = spray), data = InsectSprays) +
geom_boxplot(outlier.shape = 4) +
coord_flip()

InsectSprays[sample(nrow(InsectSprays), 10), ]
##    count spray
## 39    12     D
## 35     1     C
## 50     5     E
## 36     4     C
## 72    13     F
## 71    24     F
## 41     4     D
## 61    11     F
## 54     6     E
## 44     5     D

insects.aov <- aov(count ~ spray, data = InsectSprays)
summary(insects.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)
## spray        5   2669   533.8    34.7 <2e-16 ***
## Residuals   66   1015    15.4
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

x <- rnorm(150)
y <- rnorm(50, 0.8)
boxplot(list(x, y),
names = c("x", "y"), horizontal = TRUE,
col = 2:3, notch = TRUE, varwidth = TRUE
)
ggplot(
data = data.frame(
num = c(x, y),
idx = c(rep("x", 150), rep("y", 50))
),
aes(y = num, fill = idx)
) +
geom_boxplot(notch = TRUE) +
coord_flip()
# Wilcoxon检验的P值
wilcox.test(x, y)$p.value ## [1] 7.873804e-09 最后我们再以一个模拟数据的例子展示箱线图凹槽的功能。这里我们分别从正态分布 $$\mathrm{N}(0,1)$$$$\mathrm{N}(0.5,1)$$ 中各自产生 150 和 50 个随机数，然后作箱线图比较两组数据中间位置的差异。图 4.4 为一次模拟的结果，图中的凹槽表明了两组数据的中位数有显著差异，Wilcoxon 秩和检验也证实了这一结论。此外，该图还使用了 varwidth 参数以表明两组数据样本量的大小不同。 ## 4.4 条形图 library(RColorBrewer) # 用分类调色板 par(mfrow = c(2, 1), mar = c(3, 2.5, 0.5, 0.1)) death <- t(VADeaths)[, 5:1] barplot(death, col = brewer.pal(4, "Set1")) barplot(death, col = brewer.pal(4, "Set1"), beside = TRUE, legend.text = TRUE ) # ggplot2 reshape_VADeaths <- transform( expand.grid( sex = colnames(VADeaths), age = rownames(VADeaths) ), rates = as.vector(t(VADeaths)) ) p <- ggplot(data = reshape_VADeaths, aes(x = age, y = rates, fill = sex)) p1 <- p + geom_col(position = "stack") p2 <- p + geom_col(position = "dodge") plot_grid(p1, p2, ncol = 1) 如同前面 1.5 节中曾经提到的，条形图目前是各种统计图形中应用最广泛的，但条形图所能展示的统计量比较贫乏：它只能以矩形条的长度展示原始数值，对数据没有任何概括或推断。 R 中条形图的函数为 barplot()，用法如下： usage(barplot.default) ## ## Default S3 method: ## barplot(height, width = 1, space = NULL, names.arg = NULL, legend.text = NULL, ## beside = FALSE, horiz = FALSE, density = NULL, angle = 45, col = NULL, ## border = par("fg"), main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ## xlim = NULL, ylim = NULL, xpd = TRUE, log = "", axes = TRUE, ## axisnames = TRUE, cex.axis = par("cex.axis"), cex.names = par("cex.axis"), ## inside = TRUE, plot = TRUE, axis.lty = 0, offset = 0, add = FALSE, ## ann = !add && par("ann"), args.legend = NULL, ...) 条形图的主要参数是 height，它指定了长条的长度，这个参数可以接受一个数值向量或者一个数值矩阵作为参数，前者容易理解，后者稍有些复杂，当传入一个矩阵时，条形图针对矩阵的每一列画图，若 besideFALSE，则矩阵每一列占据一条的位置，该条由若干矩形堆砌而成，这些矩形的高度对应着矩阵的行数据，若 besideTRUE，这些矩形则并排排列而非堆砌； width 可以设置条的宽度； space 用以设置条之间的间距； names.arg 为条形图的标签，即每一条的名称； legend.text 参数在 height 为矩阵时比较有用，可以用来添加图例； horiz 用以设置条形图的方向（水平或垂直）； densityangle 等参数可以参考矩形的章节（3.4 节）； plot 为逻辑值，决定是否将条形图添加到现有图形上。 4.5 下图展示了参数 besidelegend.text 的效果。该图以 1940 年弗吉尼亚州分年龄组、分地区和分性别死亡率数据 VADeaths 为基础，展示了各组之间死亡率的差异，其中堆砌的条形图容易比较各年龄组总死亡率的大小，显然年龄越高死亡率越大，而并列的条形图容易比较组内的城乡和性别差异，一般说来，男性死亡率高于女性，农村男性死亡率低于城市男性，但女性的城乡差异没有明显规律。由于人眼对长度比比例更敏感（例如在区分城乡和性别差异时，图 4.5 的上图就不如下图直观），所以我们制图时要考虑清楚我们想展示的是数据的哪一方面，即：将最关键的信息用最能激发视觉感知的形式表现出来。 VADeaths # 弗吉尼亚州死亡数据 ## Rural Male Rural Female Urban Male Urban Female ## 50-54 11.7 8.7 15.4 8.4 ## 55-59 18.1 11.7 24.3 13.6 ## 60-64 26.9 20.3 37.0 19.3 ## 65-69 41.0 30.9 54.6 35.1 ## 70-74 66.0 54.3 71.1 50.0 ## 4.5 散点图 散点图通常用来展示两个变量之间的关系，这种关系可能是线性或非线性的。图中每一个点的横纵坐标都分别对应两个变量各自的观测值，因此散点所反映出来的趋势也就是两个变量之间的关系。 R 中散点图的函数为 plot.default()，但由于 plot() 是泛型函数（参见 B.2 小节），通常我们只需要提供两个数值型向量给 plot() 即可画散点图，或者提供一个两列的矩阵或数据框。函数 plot.default() 的用法如下： usage(plot.default) ## ## Default S3 method: ## plot(x, y = NULL, type = "p", xlim = NULL, ylim = NULL, log = "", main = NULL, ## sub = NULL, xlab = NULL, ylab = NULL, ann = par("ann"), axes = TRUE, ## frame.plot = axes, panel.first = NULL, panel.last = NULL, asp = NA, ## xgap.axis = NA, ygap.axis = NA, ...) 其中若 x 是一个两列的矩阵或数据框，则无需再提供 y ，否则 xy 都必须是数值型向量；其它参数均已在 B.2 小节中介绍。 demo("alphaDemo", package = "MSG") p <- ggplot(BinormCircle, aes(V1, V2)) + theme_void() + theme(axis.ticks = element_blank(), axis.text = element_blank(), axis.title = element_blank(), plot.background = element_rect(linetype = "solid",color = "black") ) p1 <- p + geom_point(color = rgb(1, 0, 0)) p2 <- p + geom_point(color = rgb(1, 0, 0), alpha = 0.01) plot_grid(p1, p2, ncol = 2) 4.6 展示了一个人造数据的散点图：我们设计了 2 万个样本，其中有 1 万个样本点来自于两个独立的标准正态分布，另 1 万个样本点的坐标落在半径为 0.5 的圆上，最后将这 2 万个样本拼起来并打乱顺序。该数据收录在 MSG 包中，名为 BinormCircle。虽然数据只有两个变量，但我们用普通的统计模型和数值分析几乎无法找出数据的特征，例如线性回归显示两个变量 V1V2 的回归系数非常不显著： data(BinormCircle, package = "MSG") head(BinormCircle) # 数据前 6 行 ## V1 V2 ## 1 0.889 -1.764 ## 2 0.072 -0.495 ## 3 0.123 -0.180 ## 4 -0.499 0.030 ## 5 0.252 0.432 ## 6 0.450 0.218 # 回归系数以及 P 值（不显著） coef(summary(lm(V2 ~ V1, BinormCircle))) ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.006085156 0.005258164 -1.1572776 0.2471728 ## V1 0.003988851 0.007015317 0.5685917 0.5696397 换用高阶回归的结果也类似，无论回归阶数为多少，系数均不显著，这一点从数据的构造上就可以知道（理论上两个变量的相关系数为 0）。由于样本量太大，普通的散点图上点与点之间严重重叠，所以也很难看出散点图有何异常，而使用半透明色的散点图则很容易看出，在大量的数据点中，还隐藏着一个圆圈，说明有相当一部分数据分布有特殊规律。我们在网页 https://yihui.org/en/2008/09/to-see-a-circle-in-a-pile-of-sand/ 上给出了其它五种不同的解决方案，都可以从图形的角度反映出这种规律。本章 4.21 小节也将以平滑散点图的方式再回顾这批数据。 ## 4.6 关联图 (x <- margin.table(HairEyeColor, c(1, 2))) ## Eye ## Hair Brown Blue Hazel Green ## Black 68 20 15 5 ## Brown 119 84 54 29 ## Red 26 17 14 14 ## Blond 7 94 10 16 assocplot(x) chisq.test(x)$p.value # 卡方检验 P 值
## [1] 2.325287e-25

$$$\chi^{2}=\sum_{i=1}^{r}\sum_{j=1}^{c}d_{ij}^{2};\;\; d_{ij}=\frac{f_{ij}-e_{ij}}{\sqrt{e_{ij}}} \tag{4.6}$$$

R 中关联图的函数为 assocplot()，用法如下：

usage(assocplot)
## assocplot(x, col = c("black", "red"), space = 0.3, main = NULL, xlab = NULL,
##   ylab = NULL)

4.7 是关于 HairEyeColor 数据的关联图。原始数据为一个三维数组，首先我们在性别维度上将数据汇总，得到眼睛颜色（棕蓝褐绿）和头发颜色（黑棕红金）人数的列联表。我们关心的问题是头发颜色与眼睛颜色之间是否存在关联，当然我们可以马上用函数 chisq.test() 作检验，但是检验的结果非常单一，我们只能知道零假设（独立）可否被拒绝，而图 4.7 则细致展示了数据的内部信息，例如从图中我们可以清楚看到，并非所有单元格都与期望频数有很大差异，只是少数几个单元格贡献了较大的 $$\chi^{2}$$ 值，如金发碧眼、金发棕眼等；事实上，这批数据为调查数据，眼睛颜色和头发颜色都为受访者（Delaware 大学的 592 名学生）自己填写，我们观察到金发碧眼单元格的期望频数和实际频数差异甚大，据说这背后有一则有趣的故事：由于“金发碧眼”是大家公认的美的标准，因此有些学生在填问卷时故意偏向于填写“金发碧眼”，导致“金发碧眼”的实际频数严重偏高2。从图 4.7 的代码输出中我们知道，$$\chi^{2}$$ 检验可以拒绝零假设，眼睛的颜色与头发的颜色并不独立，二者之间存在某种关联关系，然而这种关联关系是由于生物或遗传原因引起还是受访者有意隐瞒自己的信息，则需要我们进一步斟酌了。

vcd 中有一个类似的关联图函数 assoc()，但功能比本节中介绍的函数要更强大，详细介绍参见

## 4.7 条件密度图

$$$\sum_{i=1}^{k}P(Y=i|X=x)=1;\;\forall x \tag{4.7}$$$

R 中条件密度图的函数为 cdplot()，它主要是基于密度函数 density() 完成条件密度的计算 ，其用法如下：

usage(graphics:::cdplot.default)
## ## Default S3 method:
## cdplot(x, y, plot = TRUE, tol.ylab = 0.05, ylevels = NULL, bw = "nrd0",
##   n = 512, from = NULL, to = NULL, col = NULL, border = 1, main = "",
##   xlab = NULL, ylab = NULL, yaxlabels = NULL, xlim = NULL, ylim = c(0, 1), ...)
usage(graphics:::cdplot.formula)
## ## S3 method for class 'formula'
## cdplot(formula, data = list(), plot = TRUE, tol.ylab = 0.05, ylevels = NULL,
##   bw = "nrd0", n = 512, from = NULL, to = NULL, col = NULL, border = 1,
##   main = "", xlab = NULL, ylab = NULL, yaxlabels = NULL, xlim = NULL,
##   ylim = c(0, 1), ..., subset = NULL)

demo("cdplotDemo", package = "MSG")

fail <-  factor(
c(2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 1,
2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1),
levels = 2:1,
labels = c("yes", "no")
)
temperature = c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70,
70, 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81
)
fail_temperature <- data.frame(fail = fail, temperature = temperature)
ggplot(fail_temperature, aes(temperature, ..count.., fill = fail)) +
geom_density(position = "fill") +
geom_point(aes(temperature, c(0.75, 0.25)[as.integer(fail)]),
shape = 21, colour = "blue", fill = "yellow") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous("fail", expand = c(0, 0), breaks = c(0.25, 0.75),
labels = c("no", "yes"),
sec.axis = dup_axis(name = NULL,
breaks = seq(0, 1, 0.2),
labels = c(0, "", 0.4, "", 0.8, ""))
) +
scale_fill_manual(values = c("red", "lightblue")) +
theme_classic() +
theme(legend.position = "none")

## 4.8 等高图/等高线

R 中等高图的函数为 contour()，同时 grDevices 包中也提供了等高线的计算函数 contourLines()，用法分别如下：

usage(contour.default)
## ## Default S3 method:
## contour(x = seq(0, 1, length.out = nrow(z)),
##   y = seq(0, 1, length.out = ncol(z)), z, nlevels = 10,
##   levels = pretty(zlim, nlevels), labels = NULL,
##   xlim = range(x, finite = TRUE), ylim = range(y, finite = TRUE),
##   zlim = range(z, finite = TRUE), labcex = 0.6, drawlabels = TRUE,
##   method = "flattest", vfont, axes = TRUE, frame.plot = axes, col = par("fg"),
##   lty = par("lty"), lwd = par("lwd"), add = FALSE, ...)
usage(contourLines)
## contourLines(x = seq(0, 1, length.out = nrow(z)),
##   y = seq(0, 1, length.out = ncol(z)), z, nlevels = 10,
##   levels = pretty(range(z, na.rm = TRUE), nlevels))

par(mar = c(4, 4, 0.2, 0.2))
library(MSG)
data(ChinaLifeEdu)
x <- ChinaLifeEdu
library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
est <- bkde2D(x, apply(x, 2, dpik))
contour(est$x1, est$x2, est$fhat, nlevels = 15, col = "darkgreen", vfont = c("sans serif", "plain"), xlab = "预期寿命", ylab = "高学历人数" ) points(x, pch = 20) est_tidy <- data.frame( life = rep(est$x1, length(est$x2)), edu = rep(est$x2, each = length(est$x1)), z = as.vector(est$fhat)
)
levels <- pretty(range(est_tidy$z, finite = TRUE), 15) ggplot(est_tidy, aes(life, edu)) + geom_contour(aes(z = z), breaks = levels, colour = "darkgreen") + geom_point(aes(Life.Expectancy, High.Edu.NO), data = x) + labs(x ="预期寿命", y = "高学历人数") + theme_bw() 4.10 利用等高图展示了一个聚类现象。数据来源于 2005 年中国统计年鉴，数值参见 MSG 包中的 ChinaLifeEdu 数据，这里使用了其中两个变量：人口预期寿命（实际数据来自 2000 年）和高学历人口数量（定义为大专以上学历人数）。首先我们对这二维变量利用 KernSmooth 进行核密度估计，得到二维核密度值（一个矩阵），然后用两个原始变量以及这个密度值矩阵作等高图。由于密度值反映的是某个位置上数据的密集程度，图 4.10 所能揭示的现象是：中国 31 省市自治区在人口预期寿命和高学历人口数量上呈现出聚类的特征，图中密度值大的区域主要有中部、右上和左下三个，东中西格局比较明显，即：东部地区分布在图中右上角，中部省市分布在图中中部，西部地区集中在图中的左下角，对照图 4.32 可以知道聚类的具体地区名称，就更能理解这里“聚类”的含义了。关于这批数据的分析，我们在 4.26 小节仍会继续，这里不再深入。 par(mar = rep(0, 4)) # 继续前面的例子 persp(est$x1, est$x2, est$fhat, shade = 0.75, border = NA,
col = "lightblue", phi = 20, theta = 15, box = FALSE)

graphics 包中还有一个类似的等高图函数 filled.contour()，它的原理完全类似，只是它用颜色来区分高度值的大小并且有颜色图例，看起来可能更美观一些，4.12 小节中我们会详细介绍。

## 4.9 条件分割图

R 中条件分割图的函数为 coplot()，其用法如下：

usage(coplot)
## coplot(formula, data, given.values, panel = points, rows, columns,
##   show.given = TRUE, col = par("fg"), pch = par("pch"),
##   bar.bg = c(num = gray(0.8), fac = gray(0.95)),
##   xlab = c(x.name, paste("Given :", a.name)),
##   ylab = c(y.name, paste("Given :", b.name)), subscripts = FALSE,
##   axlabels = function(f) abbreviate(levels(f)), number = 6, overlap = 0.5,
##   xlim, ylim, ...)
usage(co.intervals)
## co.intervals(x, number = 6, overlap = 0.5)

co.intervals(1:10, number = 5, overlap = 0.5)
##      [,1] [,2]
## [1,]  0.5  3.5
## [2,]  2.5  5.5
## [3,]  3.5  7.5
## [4,]  5.5  8.5
## [5,]  7.5 10.5

par(mar = rep(0, 4), mgp = c(2, .5, 0))
library(maps)
coplot(lat ~ long | depth,
data = quakes, number = 4,
ylim = c(-45, -10.72), panel = function(x, y, ...) {
map("world2",
regions = c("New Zealand", "Fiji"),
add = TRUE, lwd = 0.1, fill = TRUE, col = "lightgray"
)
text(180, -13, "Fiji", adj = 1)
text(170, -35, "NZ")
points(x, y, col = rgb(0.2, 0.2, 0.2, .5))
}
)

4.12 展示了斐济岛（Fiji）附近的地震数据 quakes，数据包括地震发生地点的经纬度和震源的深度，我们想知道该地区在地震深度上分布是否均匀，因此我们令深度变量为条件变量，看在不同条件下地震发生地点（经纬度）是否有变化。从图中可以清楚看出，随着深度值的增加，地震发生地点逐渐由西向东、由南向北移动，震源较深的地震都发生在离斐济岛很近的东南侧。另外，图 4.12 还展示了 panel 参数的用法，我们借助 maps 在散点图上添加了新西兰和斐济岛的地图作为辅助信息，关于 R 中地图的使用请参考 4.34 小节。

## 4.10 一元函数曲线图

par(par(mar = c(4.5, 4, 0.2, 0.2)), mfrow = c(2, 1))
chippy <- function(x) sin(cos(x) * exp(-x / 2))
curve(chippy, -8, 7, n = 2008, xlab = "$x$", ylab = "$\\mathrm{chippy}(x)$")
curve(sin(x) / x, from = -20, to = 20, n = 200,
xlab = "$t$", ylab = "$\\varphi_{X}(t)$")

R 中函数曲线图的函数为 curve()，其用法如下：

usage(curve)
## curve(expr, from = NULL, to = NULL, n = 101, add = FALSE, type = "l",
##   xname = "x", xlab = xname, ylab = NULL, log = NULL, xlim = NULL, ...)

4.13 给出了函数 $$f(x)=\sin(\cos(x)*\exp(-x/2))$$ 的曲线以及均匀分布 $$U(-1,1)$$ 的特征函数曲线作为示例，其中特征函数为概率论中定义的 $$\varphi_{X}(t) = \mathrm{E}[\exp(itX)]$$。由于 curve() 与数据分析关系不甚密切，我们在此只是粗略介绍一下。

## 4.11 Cleveland 点图

par(mar = c(4, 4, 0.2, 0.2))
dotchart(t(VADeaths)[, 5:1], col = brewer.pal(4, "Set1"), pch = 19, cex = .65)

R 中点图的函数为 dotchart()，用法如下：

usage(dotchart)
## dotchart(x, labels = NULL, groups = NULL, gdata = NULL, offset = 1/8,
##   ann = par("ann"), xaxt = par("xaxt"), frame.plot = TRUE, log = "",
##   cex = par("cex"), pt.cex = cex, pch = 21, gpch = 21, bg = par("bg"),
##   color = par("fg"), gcolor = par("fg"), lcolor = "gray",
##   xlim = range(x[is.finite(x)]), main = NULL, xlab = NULL, ylab = NULL, ...)

4.14 再次以弗吉尼亚死亡率数据为例，给出了点图的展示。对比图 4.5 不难发现点图与条形图的相通之处。相比之下，点图的图形元素更加简洁，制图时不会显得太拥挤，我们可以视情况在这二者选其一作为表达工具。

## 4.12 颜色等高图/层次图

demo("volcano", package = "MSG")

R 中的颜色等高图函数为 filled.contour()，其用法如下：

usage(filled.contour)
## filled.contour(x = seq(0, 1, length.out = nrow(z)),
##   y = seq(0, 1, length.out = ncol(z)), z, xlim = range(x, finite = TRUE),
##   ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE),
##   levels = pretty(zlim, nlevels), nlevels = 20,
##   color.palette = function(n) hcl.colors(n, "YlOrRd", rev = TRUE),
##   col = color.palette(length(levels) - 1), plot.title, plot.axes, key.title,
##   key.axes, asp = NA, xaxs = "i", yaxs = "i", las = 1, axes = TRUE,
##   frame.plot = axes, ...)

4.15 描绘了新西兰 Maunga Whau 火山的地理数据 volcano，这份数据包含了在 $$10\mathrm{m}\times10\mathrm{m}$$ 的地理网格上测得的火山高度，是一个 $$87\times61$$ 的矩阵。仔细观察图 4.15，由于火山口的存在，颜色等高图的中部（偏左）有一小块区域的颜色并非白色，意即此处的高度比周围一圈要低。这种情况在三维图中有时未必能够迅速看出来，必须将视角调整为略向下俯视才能看到火山口。注意本图的调色板用的是“绿黄棕白”调色板，如 3.1.1 小节所介绍的，这种调色板比较适合展示地理数据。图 4.23 提供了真实的火山立体图形。

## 4.13 四瓣图

4.1 是一个典型的二维列联表，通常我们想检验的是行变量与列变量是否独立；前面 4.6 小节中曾经利用 $$\chi^{2}$$ 检验构造了关联图，这里我们从优比（Odds Ratio，OR）的角度出发对列联表进行检验。

$$$\frac{P_{1}}{P_{1}+P_{3}};\;\;\frac{P_{2}}{P_{2}+P_{4}} \tag{4.8}$$$

$$$\mathrm{OR}\equiv\frac{P_{1}}{P_{1}+P_{3}}/\frac{P_{2}}{P_{2}+P_{4}}=\frac{P_{1}(P_{2}+P_{4})}{P_{2}(P_{1}+P_{3})} \tag{4.9}$$$

$$P_{2}$$　　$$P_{4}$$ $$P_{2}+P_{4}$$
$$P_{1}+P_{2}$$　$$P_{3}+P_{4}$$

$$$\mathrm{OR}\approx\Psi\equiv\frac{P_{1}P_{4}}{P_{2}P_{3}} \tag{4.10}$$$

$$$\mathrm{Var}(\mathrm{log}(\Psi))=\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d} \tag{4.11}$$$

$$$\mathrm{log}(\Psi)\pm q_{1-\alpha/2}\sqrt{\mathrm{Var}(\mathrm{log}(\Psi))} \tag{4.12}$$$

R 中四瓣图函数 fourfoldplot() 的用法如下：

usage(fourfoldplot)
## fourfoldplot(x, color = c("#99CCFF", "#6699CC"), conf.level = 0.95,
##   std = c("margins", "ind.max", "all.max"), margin = c(1, 2), space = 0.2,
##   main = NULL, mfrow = NULL, mfcol = NULL)

ftable(UCBAdmissions) # 以紧凑形式展现出来的 UCB 录取数据
##                 Dept   A   B   C   D   E   F
## Admitted Male        512 353 120 138  53  22
##          Female       89  17 202 131  94  24
## Rejected Male        313 207 205 279 138 351
##          Female       19   8 391 244 299 317
fourfoldplot(UCBAdmissions, mfcol = c(2, 3)) # 2 行 3 列排版

4.16 是加州伯克利分校（UCB）录取数据的四瓣图，数据见于图中的代码输出。这批数据为一个 $$2\times2\times6$$ 的数组，我们可以分系别来看学生的录取是否与性别有关。从图中反映的情况来看，只有 A 系的四瓣图显示出了置信区间弧线不相交的情况，说明 A 系学生的录取与性别不独立，而其它系都不能拒绝零假设“录取与性别无关”。

(x <- apply(UCBAdmissions, 3, function(x) (x[1, 1] * x[2, 2]) / (x[1, 2] * x[2, 1])))
##         A         B         C         D         E         F
## 0.3492120 0.8025007 1.1330596 0.9212838 1.2216312 0.8278727

C 和 E 系的优比大于 1，观察图 4.16 可知， color 参数的第一个颜色值填充第一、三象限的扇形，而优比小于 1 时，第一个颜色值填充第二、四象限。

y <- qnorm(0.975) * sqrt(apply(UCBAdmissions, 3, function(x) {
sum(1 / x)
}))
conf <- exp(cbind(log(x) - y, log(x) + y))
colnames(conf) <- c("2.5%", "97.5%")
conf
##        2.5%     97.5%
## A 0.2086756 0.5843954
## B 0.3403815 1.8920166
## C 0.8545328 1.5023696
## D 0.6863345 1.2366620
## E 0.8250748 1.8087848
## F 0.4552059 1.5056335

# 对每个系的录取数据进行卡方检验分别得到 P 值
round(apply(UCBAdmissions, 3, function(x) chisq.test(x)$p.value), 3) ## A B C D E F ## 0.000 0.771 0.426 0.638 0.369 0.640 ## 4.14 颜色图 par(mar = rep(0, 4)) x <- matrix(sample(24), 8) image(1:8, 1:3, x, col = heat.colors(24), axes = FALSE, ann = FALSE) text(rep(1:8, 3), rep(1:3, each = 8), as.vector(x)) 颜色图（Color Image）与颜色等高图看起来非常类似，但是等高图需要从网格矩阵中计算等高的数据点，有时还需要一些平滑处理，而颜色图并不涉及任何背后的计算，只是简单将一个网格矩阵映射到指定的颜色序列上、以颜色方块表示数据的大小。在数据规律性较强且数据量较大的时候，这两种图形的区别可以说微乎其微，而当数据没什么规律或者数据量比较小的时候，颜色图的色块就可以很清楚地显露出来了，图 4.17 为一个简单的示意图。 R 中颜色图的函数为 image()，其用法如下： usage(image.default) ## ## Default S3 method: ## image(x = seq(0, 1, length.out = nrow(z)), y = seq(0, 1, length.out = ncol(z)), ## z, zlim = range(z[is.finite(z)]), xlim = range(x), ylim = range(y), ## col = hcl.colors(12, "YlOrRd", rev = TRUE), add = FALSE, xaxs = "i", ## yaxs = "i", xlab, ylab, breaks, oldstyle = FALSE, useRaster, ...) 参数 xyz 与等高线的参数类似，不过由于该函数为泛型函数，因此也可以接受不同类型的参数，这三个参数除了可以接受两个数值向量和一个矩阵之外， x 还可以接受一个列表，列表中包含三个子对象：x$xx$yx$z，这三个子对象分别为两个数值向量和一个矩阵，这种情况下就不需要另外单独提供 yz 参数了； col 设置一个颜色序列以便映射到不同大小的数值； add 为逻辑值，决定是否将颜色图添加到现有图形上； breaks 给定 z 分段的区间端点。

par(mar = rep(0, 4), ann = FALSE)
x <- 10 * (1:nrow(volcano))
y <- 10 * (1:ncol(volcano))
image(x, y, volcano, col = terrain.colors(100), axes = FALSE)
contour(x, y, volcano,
levels = seq(90, 200, by = 5),
add = TRUE, col = "peru"
)
box()

lattice 中提供了一个类似的函数 levelplot()，展示方法更为灵活，感兴趣的读者请参考函数的帮助文件。

## 4.15 矩阵图、矩阵点、矩阵线

sines <- outer(1:20, 1:4, function(x, y) sin(x / 20 * pi * y))
par(mar = c(2, 4, .1, .1))
matplot(sines, type = "b", pch = 21:24, col = 2:5, bg = 2:5)
# 数据矩阵的前 6 行
round(head(sines), 5)
##         [,1]    [,2]    [,3]     [,4]
## [1,] 0.15643 0.30902 0.45399  0.58779
## [2,] 0.30902 0.58779 0.80902  0.95106
## [3,] 0.45399 0.80902 0.98769  0.95106
## [4,] 0.58779 0.95106 0.95106  0.58779
## [5,] 0.70711 1.00000 0.70711  0.00000
## [6,] 0.80902 0.95106 0.30902 -0.58779

R 中矩阵图的函数为 matplot()，矩阵点的函数为 matpoints()，矩阵线的函数为 matlines()，它们的用法如下：

usage(matplot)
## matplot(x, y, type = "p", lty = 1:5, lwd = 1, lend = par("lend"), pch = NULL,
##   col = 1:6, cex = NULL, bg = NA, xlab = NULL, ylab = NULL, xlim = NULL,
##   ylim = NULL, log = "", ..., add = FALSE, verbose = getOption("verbose"))
usage(matpoints)
## matpoints(x, y, type = "p", lty = 1:5, lwd = 1, pch = NULL, col = 1:6, ...)
usage(matlines)
## matlines(x, y, type = "l", lty = 1:5, lwd = 1, pch = NULL, col = 1:6, ...)

4.19 展示了一个正弦值矩阵的矩阵图。从图上方的代码中我们可以看到矩阵的前 6 行数值，该例中只给出了参数 x 而没有 y，所以 matplot() 用矩阵 sines 的每一列依次和 1:20 画曲线。

## 4.16 马赛克图

$$$\pi_{ij}=\pi_{i\cdot}\pi_{\cdot j} \tag{4.13}$$$

$$$\mathrm{log}(\pi_{ij})=\mathrm{log}(\pi_{i\cdot})+\mathrm{log}(\pi_{\cdot j}) \tag{4.14}$$$

$$$\mathrm{log}(\mu_{ij})=\lambda+\lambda_{i}^{r}+\lambda_{j}^{c} \tag{4.15}$$$

$$\lambda_{i}^{r}$$$$\lambda_{j}^{c}$$ 分别表示行效应和列效应，$$\lambda$$ 为常数。 (4.15) 式就是最普通的对数线性模型，通过计算拟合，我们可以得到行列效应的估计值。对数线性模型在马赛克图中的主要表现是单元格的残差，而单元格的残差可以有三种：似然比残差（离差，deviance）$$G^{2}$$、Pearson $$\chi^{2}$$ 残差和 Freeman-Tukey 残差，前两种定义如下式：

$$$G^{2}=2\sum n_{ij}\mathrm{log}(\frac{n_{ij}}{\hat{\mu}_{ij}});\;\;\chi^{2}=\sum\frac{\left(n_{ij}-\hat{\mu}_{ij}\right)^{2}}{\hat{\mu}_{ij}} \tag{4.16}$$$

R 中马赛克图的函数为 mosaicplot() ，其用法如下：

usage(graphics:::mosaicplot.default)
## ## Default S3 method:
## mosaicplot(x, main = deparse1(substitute(x)), sub = NULL, xlab = NULL,
##   ylab = NULL, sort = NULL, off = NULL, dir = NULL, color = NULL,
##   shade = FALSE, margin = NULL, cex.axis = 0.66, las = par("las"),
##   border = NULL, type = c("pearson", "deviance", "FT"), ...)
usage(graphics:::mosaicplot.formula)
## ## S3 method for class 'formula'
## mosaicplot(formula, data = NULL, ..., main = deparse1(substitute(data)),
##   subset, na.action = stats::na.omit)

ftable(Titanic)
##                    Survived  No Yes
## Class Sex    Age
## 1st   Male   Child            0   5
##       Female Child            0   1
## 2nd   Male   Child            0  11
##       Female Child            0  13
## 3rd   Male   Child           35  13
##       Female Child           17  14
## Crew  Male   Child            0   0
##       Female Child            0   0
##              Adult            3  20

par(mar = c(2, 3.5, .1, .1))
mosaicplot(Titanic, shade = TRUE, main = "")

4.20 以马赛克图的形式将这个 $$4\times2\times2\times2$$ 的列联表数据展示在了同一张图中，通过矩形块（马赛克）的大小，我们可以清楚看出各舱位、不同性别、年龄的人群的生还状况。例如，对头等舱来说，无论是大人小孩或男女，下方的矩形都比上方的矩形要高（尤其是女性和小孩），这说明头等舱的生还率相对来说都比较高，很可能当时的救援是偏向头等舱的；从年龄来说，头等舱和二等舱中小孩的生存率要远高于大人，但三等舱中小孩的生存率和大人相比差异并不是太显著；但从性别角度来看，各舱位基本上还是将生存机会优先让给女性了，男性的生还率在各舱位来说都相对较低。类似地，我们还可以从图中挖掘出更多的现象，这里不再深入。另外，图中用不同颜色表示出了个单元格的残差大小，其中虚线框表示残差为负数，我们可以清楚看出哪些单元格的拟合欠佳。感兴趣的读者还可以使用 stats 包中的 loglin() 函数拟合对数线性模型、从统计模型的角度继续分析。

## 4.17 散点图矩阵

R 中散点图矩阵的函数为 pairs() ，其用法如下：

usage(pairs.default)
## ## Default S3 method:
## pairs(x, labels, panel = points, ..., horInd = 1:nc, verInd = 1:nc,
##   lower.panel = panel, upper.panel = panel, diag.panel = NULL,
##   text.panel = textPanel, label.pos = 0.5 + has.diag/3, line.main = 3,
##   cex.labels = NULL, font.labels = 1, row1attop = TRUE, gap = 1, log = "",
##   horOdd = !row1attop, verOdd = !row1attop)
usage(graphics:::pairs.formula)
## ## S3 method for class 'formula'
## pairs(formula, data = NULL, ..., subset, na.action = stats::na.pass)

4.21 是对鸢尾花数据 iris 所作的散点图矩阵，注意其中的上三角和下三角作图函数是如何定义的。对角线窗格显示的是自定义的直方图，定义如下：

# 观察如何使用 hist() 做计算并用 rect() 画图
panel.hist <- function(x, ...) {
usr <- par("usr")
on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5))
h <- hist(x, plot = FALSE)
rug(faithfuleruptions) R 中坐标轴须的函数为 rug() ，其用法如下： usage(rug) ## rug(x, ticksize = 0.03, side = 1, lwd = 0.5, col = par("fg"), ## quiet = getOption("warn") < 0, ...) 其中 x 为一个向量，给出短须的位置； ticksize 为短须的长度； side 为欲画短须的坐标轴的位置（与函数 axis() 的参数相同，参见 3.8 小节）； lwdcol 分别设定短须的宽度和颜色。 4.25 为坐标轴须的一个展示；在 4.1 小节中我们使用喷泉数据制作了直方图和密度图以便查看间隔时间的分布，这里我们在密度曲线的基础上添加上坐标轴须，以使我们能更清楚了解间隔时间数据的具体位置，而不是仅仅看一条曲线。从坐标轴须的疏密我们也可以知道分布密度的大小，这与密度曲线是相辅相成的。 严格来说，坐标轴须不能算作是图形，它只是图形的附加物，从 R 的角度来说它也只是低层作图函数，但由于它在表达数据上的优势，这里我们也将其列入本书的图库中。 ## 4.21 平滑散点图 平滑散点图的基础是散点图，但它并不直接将散点画出来，而是基于二维核密度估计 用特定颜色深浅表示某个位置的密度值大小，默认颜色越深说明二维密度值越大，即该处数据点越密集。二维核密度估计的原理和一维情况类似：一维核密度估计是在直线上按照距离远近对每个数据点加权，距离越近则对密度值贡献越大，因此数据点密集的地方的核密度值也相应大，二维情况下只是距离的计算放到了平面上，加权思想相同。 由于平滑散点图大致保留了原始数据点的位置，因此两个变量之间的关系仍然可以从图中看出来，这一点和普通的散点图类似。平滑散点图进一步的优势在于它同时还显示了二维变量的密度，从密度中我们也许可以观察到局部的聚类现象（大块的深色）。 par(mar = c(4, 4, 0.3, 0.1)) smoothScatter(BinormCircle) library(ggpointdensity) ggplot(data = BinormCircle, aes(x = V1, y = V2)) + geom_pointdensity(adjust = 0.1) 从计算的角度来说，我们首先将平面划分为 $$n\times n$$ 的网格，计算每个网格点上的二维密度值，然后用颜色将密度值大小表达出来。关于这一点，可将图 4.26 放大 8 倍，就立刻理解“网格”的意思了。 R 中平滑散点图的函数为 smoothScatter() ，其用法如下： usage(smoothScatter) ## smoothScatter(x, y = NULL, nbin = 128, bandwidth, ## colramp = colorRampPalette(c("white", blues9)), nrpoints = 100, ## ret.selection = FALSE, pch = ".", cex = 1, col = "black", ## transformation = function(x) x^0.25, postPlotHook = box, xlab = NULL, ## ylab = NULL, xlim, ylim, xaxs = par("xaxs"), yaxs = par("yaxs"), ...) 其中 xy 是两个数值向量，或者如果不提供 y 的话，可以提供一个两列的矩阵/数据框等给 xnbin 为横纵坐标方向上划分网格的数目，可以是长度为 1 或 2 的整数向量； bandwidth 为计算核密度估计时使用的带宽； colramp 为生成颜色向量的函数，默认生成从白色到蓝色渐变的颜色向量； nrpoints 为需要画出来的点的数目，因为平滑散点图的目的不是画散点，而是画颜色块，但有时候图形中某些地方的密度估计非常低，因此对应颜色也非常浅，导致读者难以察觉那些地方还有数据点的存在，此时不妨直接将这些“离群点”直接画出来，注意画点时按密度值从小到大的顺序画，一直画到第 nrpoints 个点；其它参数几乎都是画图的一般参数，此处不再介绍。 4.26MSG 包中 BinormCircle 数据的平滑散点图，其原始散点图参见图 4.6 左图。由于使用了核密度估计，这批数据中藏着的圆圈很容易就能被发现，因为它们的密度值都很大，导致这一圈上颜色也就相应较深。平滑散点图看起来和图 4.6 右图比较相似，但前者蕴含了更多的数理统计背景。不过，我们也不必一味追求数学理论，透明色可叠加这一点性质体现出的原理又何尝不是一种密度估计呢？ ## 4.22 棘状图 棘状图（Spine Plot/Spinogram）3可以看作是马赛克图（4.16 小节）的特例，也可以看作是堆砌条形图（4.4 小节）的推广。棘状图的外观看起来像是高低不齐的荆棘丛，因此而得名。 棘状图的原理和条件密度图非常类似，都展示了在给定某个自变量的情况下因变量的概率分布，但是棘状图首先对连续型的自变量进行了离散化处理，然后在离散的区间内计算因变量的条件分布。除此之外，棘状图还兼顾了自变量的分布，在横轴方向上以不同宽度的矩形表示自变量的分布密度。因此，从纵轴方向上看，棘状图用堆砌的矩形块表示因变量的分布密度，这使得它看起来像堆砌的条形图，而从横轴方向上看，棘状图用不同的矩形宽度表示自变量的密度分布，这使得它又像是马赛克图（马赛克图中矩形的长宽和频数成比例）。尤其是自变量为分类变量时，棘状图与马赛克图几乎无异。从概念和定义上来讲，棘状图实际上是用 $$P(Y|X)$$$$P(X)$$ 所作的图。 R 中棘状图的函数为 spineplot() ，其用法如下： usage(graphics:::spineplot.default) ## ## Default S3 method: ## spineplot(x, y = NULL, breaks = NULL, tol.ylab = 0.05, off = NULL, ## ylevels = NULL, col = NULL, main = "", xlab = NULL, ylab = NULL, ## xaxlabels = NULL, yaxlabels = NULL, xlim = NULL, ylim = c(0, 1), axes = TRUE, ## ...) usage(graphics:::spineplot.formula) ## ## S3 method for class 'formula' ## spineplot(formula, data = NULL, breaks = NULL, tol.ylab = 0.05, off = NULL, ## ylevels = NULL, col = NULL, main = "", xlab = NULL, ylab = NULL, ## xaxlabels = NULL, yaxlabels = NULL, xlim = NULL, ylim = c(0, 1), axes = TRUE, ## ..., subset = NULL, drop.unused.levels = FALSE) 棘状图函数是泛型函数，可以直接接受数据或者公式作为参数。 x 可以是一个分类或数值向量（自变量），也可以直接输入一个列联表，后一种情况下则不需要再输入 y 参数； y 为因变量，是一个分类变量； breaks 在自变量是连续变量时给定自变量的分段区间或者分段方法，这个参数最终被传递给直方图函数 hist() （回顾 4.1 小节）以便计算自变量的密度； off 参数指定矩形竖条之间的间距，对于连续自变量来说， off 默认为 0，对于离散自变量， off 默认为 2； col 用来设定不同类别的 y 的颜色。 # 数据取自条件密度图一节 par(mar = c(4, 4, .5, 2)) t(x <- spineplot(fail ~ temperature, col = c("lightblue", "red"))) ## temperature ## fail [50,55] (55,60] (60,65] (65,70] (70,75] (75,80] (80,85] ## yes 1 2 1 2 1 0 0 ## no 0 0 0 8 3 4 1 4.27 重新使用航天飞机 O 型环失效数据作了棘状图。我们可以看到，自变量温度被分为了 7 个长度为 5 的区间，每个区间内失效的比率都用堆砌的矩形表达。注意棘状图的横坐标轴比较特殊，它对自变量的值来说并不一定是均匀的，因为图中矩形的宽需要表达自变量的密度，所以横坐标的作用仅仅是作为密度值大小的参照物。棘状图会返回一个汇总表，表中给出了自变量的分段情况以及相应的频数。读者可以结合频数表理解棘状图。 ## 4.23 星状图、蛛网图、雷达图 # 预设调色板，stars() 默认用整数来表示颜色 palette(rainbow(12, s = 0.6, v = 0.75)) stars(mtcars[, 1:7], len = 0.8, key.loc = c(14, 1.5), ncol = 7, main = "", draw.segments = TRUE) palette("default") # 恢复默认调色板 星状图（Star Plot）、蛛网图（Spider Plot）和雷达图（Radar Plot）本质上是一类图形，它们都用线段离中心的长度来表示变量值的大小，这三种图形名称的区别在于星状图用来展示很多个多变量个体，各个个体的图形相互独立，从而整幅图形看起来就像很多星星，而蛛网图和雷达图将多个多变量个体放在同一张图形上，看起来就像是蛛网或雷达的形状，这样重叠的图形就称为蛛网图或者雷达图。简单说来，就是星状图有若干个中心，而蛛网图和雷达图只有一个中心。 R 中星状图的函数为 stars() ，其用法如下： usage(stars) ## stars(x, full = TRUE, scale = TRUE, radius = TRUE, labels = dimnames(x)[[1L]], ## locations = NULL, nrow = NULL, ncol = NULL, len = 1, key.loc = NULL, ## key.labels = dimnames(x)[[2L]], key.xpd = TRUE, xlim = NULL, ylim = NULL, ## flip.labels = NULL, draw.segments = FALSE, col.segments = 1L:n.seg, ## col.stars = NA, col.lines = NA, axes = FALSE, frame.plot = axes, main = NULL, ## sub = NULL, xlab = "", ylab = "", cex = 0.8, lwd = 0.25, lty = par("lty"), ## xpd = FALSE, ## mar = pmin(par("mar"), 1.1 + c(2 * axes + (xlab != ""), 2 * axes + (ylab != ""), 1, 0)), ## add = FALSE, plot = TRUE, ...) 参数 x 为一个多维数据矩阵或数据框，每一行数据将生成一个星形； full 为逻辑值，决定了是否使用整圆（或半圆）； scale 决定是否将数据标准化到区间 $$[0,1]$$ 内； radius 决定是否画出半径； labels 为每个个体的名称，默认为数据的行名； locations 以一个两列的矩形给出每个星形的放置位置，默认放在一个规则的矩形网格上，若提供给该参数一个长度为 2 的向量，那么所有的星形都将被放在该坐标上，从而形成蛛网图或雷达图； nrowncol 分别给定网格的行数和列数以便摆放星形，默认 nrow 等于 ncollen 为半径和线段的缩放倍数； key.loc 提供比例尺的坐标位置； key.labels 为比例尺的标签，默认为变量名称； key.xpd 设定比例尺的作图范围，参见 B.1 小节（par('xpd')）； flip.labels 设定每个星形底部的名称是否互相上下错位，以免名称太长导致文本之间互相重叠； draw.segments 设定是否作线段图，即：每个变量以一个扇形表示； col.segments 设定每个扇形区域的颜色（当 draw.segmentsFALSE 时无效）； col.stars 设定每个星形的颜色（当 draw.segmentsTRUE 时无效）； axes 决定是否画坐标轴； frame.plot 决定是否画整个图形的边框； add 决定是否将图形添加到当前图形上。 4.28 为数据 mtcars 的星状图，一共使用了 7 个变量：mpg 为每加仑汽油可行驶英里数，cyl 为汽缸数，disp 为汽缸排量，hp 为马力，drat 为背齿轮比，wt 为车重，qsec 为行驶 $$1/4$$ 英里的时间。从图中可以看到各种品牌和型号的汽车在这 7 方面的指标和性能表现。以星状图展示数据可以让我们很快找到一些有突出特征的个体，从而省去了在大批数据中逐个寻找、排序的过程。 # mtcars 数据前 7 列的前 6 行 head(mtcars[, 1:7]) ## mpg cyl disp hp drat wt qsec ## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 ## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 ## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 ## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 ## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 ## Valiant 18.1 6 225 105 2.76 3.460 20.22 ## 4.24 带状图 layout(matrix(1:2, 2), height = c(1, 1)) par(mar = c(4, 4, 0.2, 0.2)) boxplot(count ~ spray, data = InsectSprays, horizontal = TRUE, border = "red", col = "lightgreen", at = 1:6 - 0.3) stripchart(count ~ spray, data = InsectSprays, method = "stack", add = TRUE) stripchart(count ~ spray, data = InsectSprays, method = "jitter") 带状图（Strip Chart），又叫一维散点图（1-D Scatter Plot），是针对一维数据的散点图，它本质上是数据与固定值（固定 x 或固定 y）之间的散点图，这样形成的图形外观是带状的，因此称之为带状图。 R 中带状图的函数为 stripchart() ，其用法如下： usage(graphics:::stripchart.default) ## ## Default S3 method: ## stripchart(x, method = "overplot", jitter = 0.1, offset = 1/3, ## vertical = FALSE, group.names, add = FALSE, at = NULL, xlim = NULL, ## ylim = NULL, ylab = NULL, xlab = NULL, dlab = "", glab = "", log = "", ## pch = 0, col = par("fg"), cex = par("cex"), axes = TRUE, frame.plot = axes, ## ...) usage(graphics:::stripchart.formula) ## ## S3 method for class 'formula' ## stripchart(x, data = NULL, dlab = NULL, ..., subset, na.action = NULL) 带状图函数为泛型函数，可以直接接受数据参数或者公式参数。 x 为数据，一般为一个向量； method 指定作图方法，取值 overplot 意思是将所有的数据点画在一条直线上，不管它们是否有重叠，jitter 意思是将直线上的数据随机打乱，以免数据重叠导致我们不知道在某个位置究竟有多少个点，stack 意思是将重叠的数据堆砌起来，某个位置重叠的数据越多，则堆砌越高； jitter 给定打乱的程度，参见 jitter() 函数（图 6.1 有示例）； vertical 设定带状图的方向（横向或纵向）； group.names 为每一组数据的名称标签； add 决定是否将带状图添加到现有图形上； at 给定每条带子的位置。 带状图作为描述一维数据分布的工具也有其独特的优势，它的作图方法可以反映出原始数据的疏密，若原始数据有重叠，它也有相应的办法处理，最终使所有的数据点都能够被展示出来。图 4.29 重新使用了杀虫剂数据 InsectSprays，上图展示了堆砌的带状图，并且在图中同时放置了箱线图作为对比，如果只是观察箱线图，我们无从得知数据在若干位置有重复，只知道数据四分位点的位置，而带状图则可以告诉我们在哪些位置分别有多少数据点；下图为随机打乱的带状图，作图方法只是将 y 方向上的固定数值添加了随机数，使原本重叠在一起的数据重新拥有不同的纵坐标，从而将重叠的数据分开来。 回顾前面 4.20 小节中我们曾经用坐标轴须在坐标轴上标记出原始数据，这与带状图在一条直线上用点表达原始数据的想法有异曲同工之妙，然而坐标轴须本身没有堆砌和随机打乱数据的功能（可以通过函数 jitter() 函数实现），所以有时候使用不妨视情况向图中添加带状图作为变量分布的辅助性描述。 ## 4.25 向日葵散点图 par(mar = c(4, 4, 0.2, 0.2)) sunflowerplot(iris[, 3:4], col = "gold", seg.col = "gold") 向日葵散点图（Sunflower Scatter Plot）是用来克服散点图中数据点重叠问题的特殊散点图工具。它采用的办法是在有重叠的地方用一朵“向日葵花”的花瓣数目来表示重叠数据的个数，这样我们就很容易看出来散点图中哪些地方的数据有重叠，而且能知道重叠的具体数目。向日葵散点图在数据特别密集或者数据类型为分类数据时很有用，因为这两种情况下都容易产生重复的数据点，尤其是后一种情况下，数据几乎必然有重复（除非列联表单元格频数为 1）。 R 中向日葵散点图的函数为 sunflowerplot() ，其用法如下： usage(graphics:::sunflowerplot.default) ## ## Default S3 method: ## sunflowerplot(x, y = NULL, number, log = "", digits = 6L, xlab = NULL, ## ylab = NULL, xlim = NULL, ylim = NULL, add = FALSE, rotate = FALSE, pch = 16, ## cex = 0.8, cex.fact = 1.5, col = par("col"), bg = NA, size = 1/8, ## seg.col = 2, seg.lwd = 1.5, ...) xy 分别为散点图的两个变量； number 为人工给定的数据频数，即图中的花瓣数目，若不指定这个参数的话 R 会自动从 xy 计算； digits 给定数值的有效数字位数，在计算重复数据之前原始数据会按照 digits 四舍五入； add 决定是否将向日葵散点图添加到当前图形上； rotate 决定是否随机旋转向日葵的角度； pch 给定散点图的点的类型； cex 给定散点图的点的缩放倍数； cex.fact 给定向日葵中心点的缩小倍数，真正的缩放倍数为 cex/cex.factcol 为散点的颜色， bg 为点的背景色； size 为向日葵花瓣的长度，单位为英寸； seg.col 为花瓣颜色； seg.lwd 为花瓣宽度。 4.30 为鸢尾花花瓣长和宽的向日葵散点图，注意左下方和中部都有一些重复数据，这幅图的颜色用的是金色 'gold'，目的是为了使向日葵散点图看起来与其名称相符，不过似乎使用鸢尾花的颜色更符合数据的背景，读者不妨试试看。 ## 4.26 符号图 符号图是用各种符号展示高维数据的图示工具，它的主要思想是将高维数值体现在图形中符号的特征上。因为受到平面的限制，我们对于高维数据的展示方法总是很有限，仅仅是对于二维数据的展示最为方便，对更高维度的如三维、四维甚至五维的数据相对缺乏好的工具。由于符号的存在，使得我们可以将高于二维的数据“寄托”在符号的各种特征上，如：以矩形为散点图的基本符号，那么我们可以用其长宽分别代表两个变量，这样一幅图形中至少可以放置四个变量；类似地，我们可以以圆圈、正方形、多边形、箱线图、温度计等符号为散点图中的“点”，于是散点图就可以被扩展为高维数据的展示工具。 R 中的符号图函数为 symbols()，它提供了六种基本符号：圆、正方形、长方形、星形、温度计和箱线图，分别由相应的参数指定； symbols() 的用法如下： usage(symbols) ## symbols(x, y = NULL, circles, squares, rectangles, stars, thermometers, ## boxplots, inches = TRUE, add = FALSE, fg = par("col"), bg = NA, xlab = NULL, ## ylab = NULL, main = NULL, xlim = NULL, ylim = NULL, ...) 如前所述，符号图的基础是散点图，因此首先要给出两个参数 xy 以便作散点图，然后在散点的位置上画出符号；接下来的六个参数分别指定符号的形状： circles 圆：一个数值向量，给定圆的半径（实际上是与圆的半径成比例，下同） squares 正方形：一个数值向量，给定正方形的边长 rectangles 长方形：一个矩阵，列数为 2，这两列分别给定长方形的宽和高 stars 星形：一个矩阵，列数 $$\geq3$$，类似雷达图，给定从星星中心向每个方向的射线的长度（严格说是线段），最终这些射线的端点会连接起来形成一个星形，但射线本身不会被画出来；缺失值将被视作 0 （星形在符号图中并不直观，推荐直接使用星状图，见 4.28 小节） thermometers 温度计：一个矩阵，列数为 3 或 4，前两列分别给定温度计的宽和高；若矩阵为三列，那么第三列为温度计内的“温度”高度，注意这一列的值应该小于 1，否则温度的填充会超过温度计的范围；若矩阵为四列，那么温度将按照第三列与第四列的比率进行填充，同样，这两列的比率需要小于 1 boxplots 箱线图：一个矩阵，列数为 5，前两列分别给定箱子的宽和高，第三、四列分别给定两条线（下线和上线）的长度，第五列与温度计类似，给定箱线图内的中位数标记线在箱子内部的高度比例，因此这一列数据也需要在 $$[0,1]$$ 范围内；这里只是借用了箱线图的称谓，符号图中的箱线与真正的箱线图之间没有关系 demo("symbols_all", echo = FALSE, package = "MSG") 不难看出，这六种符号图能展示的数据维度分别为 3、3、4、$$\geq5$$、5 或 6、7；参数 inches 为逻辑值，控制着符号的大小，若为 TRUE（默认），那么图中所有符号的最长长度（边长或半径等等）将被设置为 1 英寸（约 2.54 厘米），其它长度按比例缩放，若该参数取一个正数，那么类似地，所有符号的最长长度的英寸值为该正数，因此 TRUE 和 1 的效果是相同的，如果 inchesFALSE 那么符号的长度单位取相应坐标轴的单位，例如符号中的 1 长度即为坐标轴上的 1 单位； add 为逻辑值，设定是否将符号图添加到现有图形上； fg 为符号的前景色； bg 为符号的背景色或填充色；其它参数用来添加标题、设定坐标轴范围等。 4.31 给出了六种基本符号的形状，注意观察各种符号是如何利用自己的特征将高维数据表达出来的。下面我们通过实际数据来看符号图在展示数据时的效果，以 2005 年中国 31 省市自治区的人口特征数据为例，我们选取了人口自然增长率、年底人口总数、城镇人口比重、人口预期寿命和高学历人数五个变量作为国民素质的刻画指标，数据见 MSG 包中的 ChinaPop.R 文件（可用 source() 函数载入，见以下代码），这些数据都可以从《中国统计年鉴》中获得。 source(system.file("extdata", "ChinaPop.R", package = "MSG"), encoding = "UTF-8") head(ChinaPop) ## 增长率 总人口 城镇人口比重 预期寿命 高学历人数 ## 北京 1.09 1538 0.8362 76.10 48001 ## 天津 1.43 1043 0.7511 74.91 18601 ## 河北 6.09 6851 0.3769 72.54 40036 ## 山西 6.02 3355 0.4211 71.65 23197 ## 内蒙古 4.62 2386 0.4720 69.87 23660 ## 辽宁 0.97 4221 0.5870 73.34 44404 demo("ChinaPop", package = "MSG") ## Warning in (function (z, notch = FALSE, width = NULL, varwidth = FALSE, : some ## notches went outside hinges ('box'): maybe set notch=FALSE 4.32 融合了多种图形和图形元素，它的基础是一幅等高图（回顾 4.8 小节和图 4.10），利用人口预期寿命和高学历人数两个变量计算二维密度，画出等高线，便完成了底图的制作；然后我们通过人口预期寿命和高学历人数两个变量的数值往图中添加温度计符号，温度计宽代表增长率，高代表总人口数，温度代表城镇人口比重；然后我们用 text() 函数将各省市的文本标签添加到图中。经过这些图形元素的表达，全国 31 省市自治区的五项人口特征便一目了然，例如通过温度计的高度可以观察出三个人口大省广东、山东、河南（相应的人口总量小的地区如西藏、青海、宁夏等也容易看出），由宽度可以看出西藏、青海、宁夏、新疆等省市自治区的人口自然增长率非常高（而北京、上海、天津等直辖市的增长率则很低），从温度指示的情况来看，北京、上海和天津三大直辖市的城镇人口比例要远高于其它地区；从整幅散点图来看，人口平均预期寿命与高学历者人数呈比较明显的正相关关系。箱线图和坐标轴须分别刻画了人口平均预期寿命与高学历者人数各自的分布特征。这样，我们就完成了在平面上描述五维变量的任务。从这幅图我们可以看出掌握基本图形元素使用的用处和重要性 — 没有一种统计软件能够提供现成的模块或函数来完成类似的任务。 ## 4.27 饼图 layout(matrix(c(1, 2, 1, 3), 2)) # 拆分作图区域 par(mar = c(4, 4, 0.2, 0.2)) pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12) names(pie.sales) <- c("Blueberry", "Cherry", "Apple", "Boston Cream", "Other", "Vanilla Cream") pie.col <- c("purple", "violetred1", "green3", "cornflowerblue", "cyan", "white") pie.sales <- sort(pie.sales, decreasing = TRUE) # 排序有助于可读性 pie(pie.sales, col = pie.col) dotchart(pie.sales, xlim = c(0, 0.3)) barplot(pie.sales, col = pie.col, horiz = TRUE, names.arg = "", space = 0.5) 饼图是目前应用非常广泛的统计图形，然而，根据统计学家（主要是 Cleveland 和 McGill）和一些心理学家的调查结果 ，这种以比例展示数据的统计图形实际上是很糟糕的可视化方式，因此，R 关于饼图的帮助文件中清楚地说明了并不推荐使用饼图，而是使用条形图（4.4 小节）或点图（4.11 小节）作为替代。读者若有兴趣可以到 COS 论坛查看这则关于饼图的幽默：https://d.cosx.org/d/101157 饼图的原理很简单，每一个扇形的角度与相应数据的数值大小成比例，关于饼图的知识此处不再赘述。R 提供了函数 pie() 制作饼图，用法如下： usage(pie) ## pie(x, labels = names(x), edges = 200, radius = 0.8, clockwise = FALSE, ## init.angle = if (clockwise) 90 else 0, density = NULL, angle = 45, ## col = NULL, border = NULL, lty = NULL, main = NULL, ...) 参数 x 为一个数值向量， labels 为标签，其它参数基本上都是为多边形准备的，因为扇形实际上是多边形所作，例如 edges 可以设定圆弧的光滑程度（多边形的边越多则越光滑）， densityangle 等参数参见多边形的章节（3.4 小节）。 4.33 同时给出了饼图、点图和条形图分别对一个不同牌子馅饼的销售数据的展示，请读者对比并思考饼图对于展示数据的弱势，例如从饼图中看 Boston Cream 和 Vanilla Cream 差异有多大，而在条形图中看差异又是多大？两种情况下我们对差异的感受相同吗？ 至此，我们已经全部介绍完 graphics 包中的统计图形函数，读者一般不用安装别的附加包就可以完成以上的图形制作。下面我们开始选择性介绍其它基础包和附加包中的图形函数和图形种类。 ## 4.28 热图 热图的基础就是 4.14 小节中介绍的颜色图，但它在颜色图上做了一个特殊处理，就是聚类。具体来说，热图也是将一个矩阵中单元格数值用颜色表达，如颜色深表示数值大，但热图并非只是简单表达数值大小，而是对矩阵的行或列进行层次聚类，获得聚类的结果之后将行或列以聚类的顺序排列，并在颜色图的边界区域加上聚类的谱系图。这样一来，我们不仅可以直接观察矩阵中的数值分布状况，也可以立即知道聚类的结果，可谓一举两得。关于聚类分析的进一步介绍，参见 6.2.11 小节。 R 中热图函数为 stats 包中的 heatmap() ，其用法如下： usage(heatmap) ## heatmap(x, Rowv = NULL, Colv = if (symm) "Rowv" else NULL, distfun = dist, ## hclustfun = hclust, reorderfun = function(d, w) reorder(d, w), add.expr, ## symm = FALSE, revC = identical(Colv, "Rowv"), ## scale = c("row", "column", "none"), na.rm = TRUE, margins = c(5, 5), ## ColSideColors, RowSideColors, cexRow = 0.2 + 1/log10(nr), ## cexCol = 0.2 + 1/log10(nc), labRow = NULL, labCol = NULL, main = NULL, ## xlab = NULL, ylab = NULL, keep.dendro = FALSE, ## verbose = getOption("verbose"), ...) 其中 x 是数据矩阵，它的类型只能是矩阵，不能是数据框或其它类型； RowvColv 分别决定了行和列如何计算层次聚类和重新排序，若设置为 NULL（默认），则按层次聚类的结果将行和列重新排序并相应画谱系图，若为 NA 的话，则不画谱系图； distfun 决定用哪个函数计算距离以便进一步计算聚类，默认为 dist()hclustfun 决定用哪个函数计算层次聚类； scale 决定对行或列进行标准化，或者不进行标准化； margins 设定图形的下边距和右边距； ... 参数传递给 image() ，所以我们还可以利用 image() 的参数来调整图形外观，比如用 col 设置单元格的颜色系列。 # 用极端化调色板 library(RColorBrewer) heatmap(as.matrix(mtcars), col = brewer.pal(9, "RdYlBu"), scale = "column", margins = c(4, 8)) 4.34 展示了 Motor Trend 杂志 1974 年汽车数据的热图。在图 4.28 中我们曾经部分使用过这批数据并说明了变量的含义，读者可以对比该图，看用颜色和用星形哪种方式表达数值大小更易感知。在汽车数据热图 4.34 中，我们使用了极端化调色板，用以强调极端值（回顾 3.1.1 小节），所以很容易观察到各项较大或较小的汽车性能指标，如马力最大的是 Maserati Bora。从行的聚类来说，可以看到同一品牌的不同型号容易聚在一起，如 Mazda 的两款车，通过颜色的比较，我们又可以看出聚在一类的车中，差异在哪个或哪些指标上；从列的聚类来说，马力 hp 和气缸排量 disp 两个变量比较相似，聚为了一类，而且它们最后才和其它指标聚成一类，说明这两个指标和其它指标的差异较大，可以想象，如果做主成分分析，这两个指标也许可以提取一个成分。注意图中的数据是对列进行过标准化的，如果不做标准化，那么聚类结果就容易被数量级大的变量主导，导致产生一些误解。 普通颜色图只是按照原始行列顺序排列色块，所以看起来可能显得比较混乱，而通常来说热图中的色块看起来会稍微整齐一些，颜色相近的色块往往会排在相近的地方，这样一来，如果样本数据中真的存在很明显的聚类现象，那么在热图中的直接反映就是不同颜色色块的高度集中。 ## 4.29 交互效应图 在回归模型或方差分析中，我们常遇到交互效应的概念。所谓交互效应，就是一个自变量对因变量的影响大小受另一个变量取值水平的影响，以二元回归为例，以下就是一个典型的含有交互效应的回归模型： $y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{12}x_{1}x_{2}+\epsilon$ $$c_{2}\equiv\beta_{2}+\beta_{12}x_{1}$$，上式可以改写为： \begin{aligned} y & = \beta_{0}+\beta_{1}x_{1}+(\beta_{2}+\beta_{12}x_{1})x_{2}+\epsilon\\ & \equiv \beta_{0}+\beta_{1}x_{1}+c_{2}x_{2}+\epsilon\end{aligned} 显然，变量 $$x_{2}$$ 对因变量 $$y$$ 的影响系数 $$c_{2}$$$$x_{1}$$ 取值水平的控制，同理 $$x_{1}$$ 对 $$y$$ 的影响也受 $$x_{2}$$ 取值水平的影响。若模型中不存在交互项 $$x_{1}x_{2}$$，那么无论 $$x_{1}$$ 取值水平如何，$$x_{2}$$ 每增加 1 单位，$$y$$ 都会变化 $$\beta_{2}$$ 单位，对 $$x_{2}$$ 而言同理。以上数学原理用更直白的话来说就是，若因变量随一个变量的变化幅度与另一个变量的取值水平无关，那么模型中不存在交互效应。这个现象反映到图中，便有了交互效应图，它通常是针对分类变量而言的，看一个分类变量给定分类水平时，因变量在另一个分类变量各水平下的均值如何变化，这种变化趋势如果在前一个分类变量换一个取值水平后仍然保持相同的话，则说明这两个分类变量没有交互效应。 R 中交互效应图的函数为 interaction.plot()，用法如下： usage(interaction.plot) ## interaction.plot(x.factor, trace.factor, response, fun = mean, ## type = c("l", "p", "b", "o", "c"), legend = TRUE, ## trace.label = deparse1(substitute(trace.factor)), fixed = FALSE, ## xlab = deparse1(substitute(x.factor)), ylab = ylabel, ## ylim = range(cells, na.rm = TRUE), lty = nc:1, col = 1, ## pch = c(1L:9, 0, letters), xpd = NULL, leg.bg = par("bg"), leg.bty = "n", ## xtick = FALSE, xaxt = par("xaxt"), axes = TRUE, ...) 其中 x.factor 是横坐标上的分类变量； trace.factor 是第二个分类变量，按照这个分类变量的不同取值水平将 x.factor 分类下的因变量均值连接起来； response 是因变量； fun 是指定的对因变量汇总的函数，默认为均值，当然我们也可以指定其它计算函数如中位数 median()type 为画图类型，见 B.2 小节图 B.4legend 决定是否画图例；其它参数用以设置图形细节，如各条均值连线的线条样式等。 par(mar = c(4, 4, 0.2, 0.2)) with(esoph, { interaction.plot(agegp, alcgp, ncases / (ncases + ncontrols), trace.label = "alcohol", fixed = TRUE, xlab = "Age", ylab = "Cancer Proportion") # 方差分析，交互项系数不显著 summary(aov(ncases / (ncases + ncontrols) ~ agegp * alcgp)) }) ## Df Sum Sq Mean Sq F value Pr(>F) ## agegp 5 3.505 0.7010 14.954 1.04e-09 *** ## alcgp 3 3.474 1.1580 24.704 9.79e-11 *** ## agegp:alcgp 15 1.130 0.0753 1.607 0.0968 . ## Residuals 64 3.000 0.0469 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 4.35 展示了一批法国食道癌数据的交互效应图，这个数据名为 esoph，在 datasets 包中，它记录了一组人群的年龄、饮酒量、抽烟量以及是否患食道癌，前三个变量都是以分组形式记录的分类变量。我们不难想象，年龄越高则患食道癌的可能性就会越大，但我们现在关心的问题是，在不同饮酒量水平下，年龄对食道癌的影响是一致的吗？图 4.35 对这个问题的回答大致是肯定的，因为图中各条折线的增长趋势是差不多的，也就是说，不管人群属于哪种饮酒量水平，年龄对患食道癌概率的影响关系都大致一样，但其中也有一些微小的差异，例如年龄组从 35-44 到 45-54，患食道癌的概率在饮酒量大于 40 克/天的人群中会增长大约 20 %，但饮酒量小于 40 克/天的人群中，患食道癌的概率反而略微有所下降，在这里我们可以看到一点交互总用，但总体来说，各饮酒组中，年龄对食道癌的影响方向还是大致一致的。这种似有似无的交互作用可以用方差分析来检验，结果如图 4.35 的代码输出，交互项的 P 值大于 0.05， 说明交互项基本上可以忽略（系数不显著）。总之，从交互效应图看两个分类变量是否有交互作用，只需要看各条折线是否大致平行即可。 我们也可以计算一下每种分组组合下的因变量均值，辅助理解图 4.35 tbl <- aggregate(ncases / (ncases + ncontrols) ~ agegp * alcgp, data = esoph, mean) colnames(tbl) <- c("agegp", "alcgp", "mean") reshape(tbl, timevar = "alcgp", idvar = "agegp", direction = "wide", sep = "_") ## agegp mean_0-39g/day mean_40-79 mean_80-119 mean_120+ ## 1 25-34 0.000000000 0.00000000 0.0000000 0.2500000 ## 2 35-44 0.017857143 0.05046584 0.0000000 0.3888889 ## 3 45-54 0.005434783 0.34899749 0.3290179 0.8541667 ## 4 55-64 0.273461657 0.31150210 0.6333333 0.7142857 ## 5 65-74 0.168898810 0.45185185 0.6153846 0.8125000 ## 6 75+ 0.240740741 0.43333333 1.0000000 1.0000000 ## 4.30 QQ 图 关于统计分布的检验有很多种，例如 KS 检验、卡方检验等，从图形的角度来说，我们也可以用 QQ 图（Quantile-Quantile Plots）来检查数据是否服从某种分布。QQ 图的原理并不复杂：如果一批数据 $$x_{1},x_{2},\ldots,x_{n}$$ 服从某种理论分布，那么将排序后的数据 $$x_{(1)},x_{(2)},\ldots,x_{(n)}$$ 和理论分布的分位数 $$q_{1/n},q_{2/n},\ldots,q{}_{n/n}$$ 去画散点图，得到的 $$n$$ 个点应该大致排列在对角线上，因为这两批数字应该大致相等。从另一个角度来看，检验一批数据是否服从某种理论分布，也就是看其经验分布和理论分布是否一致，而排序后的数据 $$x_{(1)},x_{(2)},\ldots,x_{(n)}$$ 可以看作是经验分布的 $$1/n,\,2/n,\,\cdots,n/n$$ 分位数，若这些分位数和理论分位数一致，也就说明了经验分布和理论分布相似。为了说明这一点，我们可以看看数值模拟的结果： # 从 N(0, 1) 中生成 1000 个随机数的分位数 quantile(rnorm(1000), probs = seq(.1, .9, .2)) ## 10% 30% 50% 70% 90% ## -1.26732829 -0.47817587 0.01519423 0.51877480 1.29469538 # 真实的分位数 qnorm(seq(0.1, .9, .2)) ## [1] -1.2815516 -0.5244005 0.0000000 0.5244005 1.2815516 以上数据的 5 个分位数和理论分位数都比较接近，读者可以模拟其它分布，例如从卡方分布中生成随机数，看其分位数是多少，与正态分布分位数差异如何。 R 中 QQ 图的函数为 qqplot() ，由于正态分布是我们经常检验的分布，R 也直接提供了一个画正态分布 QQ 图的函数 qqnorm() ，这两个函数都在基础包 stats 包中，它们的用法如下： usage(qqplot) ## qqplot(x, y, plot.it = TRUE, xlab = deparse1(substitute(x)), ## ylab = deparse1(substitute(y)), ...) usage(stats:::qqnorm.default) ## ## Default S3 method: ## qqnorm(y, ylim, main = "Normal Q-Q Plot", xlab = "Theoretical Quantiles", ## ylab = "Sample Quantiles", plot.it = TRUE, datax = FALSE, ...) usage(qqline) ## qqline(y, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75), ## qtype = 7, ...) 由于 qqplot() 检验的是两批数据的分布是否相同，所以它需要两个数据参数 x 和 y ， qqnorm() 只需要一个数据参数 x ，其它设置标签和标题等元素的图形参数此处不再赘述。 par(mfrow = c(1, 2)) par(mar = c(4, 4, 0.2, 0.2)) x <- scale(geyserwaiting)
qqnorm(x, cex = 0.7, asp = 1, main = "")
abline(0, 1)
plot(density(x), main = "", xlim = range(x))
curve(dnorm, from = -3, to = 3, lty = 2, add = TRUE)

4.36 左图是喷泉间隔时间数据的正态分布 QQ 图（4.1 小节的直方图用到过），注意其中的数据经过了标准化，使之均值为 0，方差为 1。可以看出，数据点并不呈直线分布，这说明（标准化后的）数据的分布和标准正态分布有所差异，那么具体是何种差异呢？图的左边有一部分点偏离在直线上方，说明实际分位数大于理论分位数，从密度曲线的角度来说，也就是实际数据的分布曲线更偏右一些，理论分布曲线左边的尾巴向左伸得更远，而图的右边又有一些点在直线下方，说明此处实际分布曲线偏左，即实际分位数偏小。右图画出了数据的核密度估计曲线（实线）和真正的标准正态分布密度曲线（虚线），读者可以将二者的对比结合左图来理解 QQ 图中数据点偏离直线的方向与分布曲线的偏向关系。

QQ 图的用途不仅在检查数据是否服从某种特定理论分布，它也可以推广到检查数据是否来自某个位置参数分布族。例如，若数据来自正态分布 $$N(5,1)$$，我们拿它和标准正态分布 $$N(0,1)$$ 画 QQ 图的话，数据点仍然会大致排列在直线上，此时直线的斜率仍然是 1，但截距就不是 0 了，而是 5。正态分布是位置参数分布族中的一种分布：若 $$X\sim N(\mu,\sigma^{2})$$，那么 $$X + \delta \sim \mathrm{N}(\mu+\delta,\sigma^{2})$$ 仍然是正态分布。均值 $$\mu$$ 是位置参数。

## 4.31 生存函数图

library(survival)
leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml)
plot(leukemia.surv, lty = 1:2, xlab = "time")
legend("topright", c("Maintenance", "No Maintenance"), lty = 1:2, bty = "n")
survminer::ggsurvplot(leukemia.surv, data = aml)

$\hat{S}(t)=\left\{ \begin{array}{ll} \prod_{i:\, t_{(i)}\leq t}(\frac{n-i}{n-i+1})^{\delta_{(i)}}, & \mbox{对}t\leq t_{(n)};\\ \left\{ \begin{array}{ll} 0 & \mbox{如果}\delta_{(n)}=1,\\ \mbox{未定义} & \mbox{如果}\delta_{(n)}=0, \end{array}\right. & \mbox{对}t>t_{(n)}. \end{array}\right.$

survival 提供了生存函数的计算和估计方法。具体函数为 survfit()，它返回一个 survfit 类的对象，而 survival 包扩展了泛型函数 plot()，使其拥有子函数 plot.survfit()，因此在估计完生存函数之后，我们可以直接调用 plot() 生成生存函数图。函数 plot.survfit() 的用法如下：

usage(survival:::plot.survfit)
## ## S3 method for class 'survfit'
## plot(x, conf.int, mark.time = FALSE, pch = 3, col = 1, lty = 1, lwd = 1,
##   cex = 1, log = FALSE, xscale = 1, yscale = 1, xlim, ylim, xmax, fun,
##   xlab = "", ylab = "", xaxs = "r", conf.times, conf.cap = 0.005,
##   conf.offset = 0.012,
##   conf.type = c("log", "log-log", "plain", "logit", "arcsin"), mark,
##   noplot = "(s0)", cumhaz = FALSE, firstx, ymin, ...)

4.37 展示了急性髓细胞白血病（Acute Myelogenous Leukemia）数据 aml 的生存函数图，该数据有一个分组变量 x 表示病人是否接受了化疗，从图中可以看出，接受化疗的病人生存函数的下降速度比没接受化疗的病人要慢，表明化疗还是有一定作用的。进一步我们可以用对数秩检验知道这两组病人的生存时间在 10% 的显著水平下有显著差异：

survdiff(Surv(time, status) ~ x, data = aml)
## Call:
## survdiff(formula = Surv(time, status) ~ x, data = aml)
##
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained    11        7    10.69      1.27       3.4
## x=Nonmaintained 12       11     7.31      1.86       3.4
##
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.07

## 4.32 分类与回归树图

rpart 提供了分类与回归树的计算拟合函数 rpart() ，该函数包同时也扩充了泛型函数 plot() ，凡是 rpart 类型的对象在作图时都会自动调用 plot.rpart() 生成树图。plot.rpart() 的用法如下：

library(rpart)
usage(rpart:::plot.rpart)
## ## S3 method for class 'rpart'
## plot(x, uniform = FALSE, branch = 1, compress = FALSE, nspace, margin = 0,
##   minbranch = 0.3, ...)

x 是一个 rpart 类型的对象，一般由 rpart() 函数拟合产生； uniform 决定是否在从上至下的枝节点之间使用相等的纵向距离以避免树枝在某些局部区域靠得太近使图形难以辨认，默认情况下每两个枝节点之间的距离与拟合误差成比例； branch 设定树枝的形状，0 为 “V” 字型，1 为垂直的形状，该参数可以取 $$[0,1]$$ 之间的数值以使得数值形状更像 “V” 或更垂直； compress 设定是否在横向上压缩树枝的间距使得图形更紧凑。

fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis)
par(mar = rep(1, 4), xpd = TRUE)
plot(fit, branch = 0.7)
text(fit, use.n = TRUE, digits = 7)

## 4.33 小提琴图

library(sm)
library(vioplot)
usage(vioplot)
## vioplot(x, ...)
f <- function(mu1, mu2)
c(rnorm(300, mu1, 0.5), rnorm(200, mu2, 0.5))
x1 <- f(0, 2)
x2 <- f(2, 3.5)
x3 <- f(0.5, 2)
vioplot(x1, x2, x3,
horizontal = TRUE, col = "bisque",
names = c("A", "B", "C")
)

4.39 用三个随机数序列展示了小提琴图的外观及其在表达数据密度和比较统计分布参数（中位数、众数等）上的功效。

lattice 包中的函数 panel.violin() 也提供了类似的小提琴图展示。关于小提琴图的理论请参考

## 4.34 地图

maps 包中核心的函数为 map() ，它的用法如下：

library(maps)
usage(map)
## map(database = "world", regions = ".", exact = FALSE, boundary = TRUE,
##   interior = TRUE, projection = "", parameters = NULL, orientation = NULL,
##   fill = FALSE, col = 1, plot = TRUE, add = FALSE, namesonly = FALSE,
##   xlim = NULL, ylim = NULL, wrap = FALSE, resolution = if (plot) 1 else 0,
##   type = "l", bg = par("bg"), mar = c(4.1, 4.1, par("mar")[3], 0.1),
##   myborder = 0.01, namefield = "name", lforce = "n", ...)

demo("AgriComp", package = "MSG")

4.40 展示了 2005 年世界各国家地区的农业进出口竞争力指标数据 (Xie 2007)，其中我们将竞争力指标简单定义为（出口-进口）/（出口+进口）。地图上方我们自行添加了颜色图例，从图中可以看出，阿根廷、巴西等南美国家的农业进出口竞争力较强，而利比亚、阿尔及利亚等北非国家的竞争力较弱。该地图的大致制作过程为：首先我们用 world 数据库作出一幅空白的世界地图，地区边界用灰色线条表示，然后我们根据竞争力数据中的地区名称与地理数据库中地区名称的对应将数据以颜色的形式表示到世界地图中，最后我们在图中添加了赤道线以及东盟国家（ASEAN）的矩形区域，因为会议论文 Xie (2007) 的主题是中澳自由贸易区。

## 4.35 脸谱图

TeachingDemos 提供了两个脸谱图函数 faces()faces2() ，两个函数能反映的面部特征不尽相同，各有所长，例如 faces() 可以画头发和耳朵，但 faces2() 可以画更多的变量，这里我们只介绍后者，读者可以阅读前者的帮助文档了解更多信息。faces2() 的用法如下：

library(TeachingDemos)
usage(faces2)
## faces2(mat, which = 1:ncol(mat), labels = rownames(mat),
##   nrows = ceiling(nrow(mat)/ncols), ncols = ceiling(sqrt(nrow(mat))),
##   byrow = TRUE, scale = c("columns", "all", "center", "none"),
##   fill = c(0.5, 0.5, 1, 0.5, 0.5, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 0.5),
##   ...)

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

faces2(mtcars[, c("hp", "disp", "mpg", "qsec", "wt")], which = c(14, 9, 11, 6, 5))

4.41 是汽车数据 mtcars 中 5 个变量的脸谱图，这 5 个变量分别为马力 hp、气缸排量 disp、每加仑行驶英里数 mpg、行驶 1/4 英里时间 qsec 和车重 wt。我们将 which 参数设定为 c(14, 9, 11, 6, 5)，也就是指定这几个变量分别用眼睛大小、嘴宽、眼距、鼻长和下半脸宽来表示，所以整幅图形的解读方式就是：眼睛瞪得越大，说明该车型的马力越大；嘴越宽则气缸排量越大；两眼距离越大则越省油（每加仑汽油跑得越远）；鼻子越长则说明跑得越快；脸越宽则说明车越重。那么，我们很容易看出，Maserati Bora 马力最强（大眼睛），Lincoln Continental 等车气缸排量较大（宽嘴），Honda Civic 等车比较省油（眼距大），Merc 230 跑得最快（鼻子长），Lincoln Continental 等车最重（脸胖）。这里我们尽量将这几个汽车性能指标形象化到合适的脸部特征上（如瞪着大眼睛表示马力大），读者在遇到具体的案例数据时，不妨也仔细考虑一下指标的实际意义以及安排它们到哪个面部表情上。很多人都知道宏基施振荣的“微笑曲线”并且对之印象深刻，其原因何尝不是因为这条曲线形象而且直观呢？

## 4.36 平行坐标图

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from
##   +.gg   ggplot2
usage(ggparcoord)
## ggparcoord(data, columns = 1:ncol(data), groupColumn = NULL, scale = "std",
##   scaleSummary = "mean", centerObsID = 1, missing = "exclude", order = columns,
##   showPoints = FALSE, splineFactor = FALSE, alphaLines = 1, boxplot = FALSE,
##   shadeBox = NULL, mapping = NULL, title = "")

ggparcoord(iris, columns = 1:4, groupColumn = 5, scale = "uniminmax") +
geom_line(size = 1.2)

4.42 是鸢尾花数据的平行坐标图，数据经过了默认的标准化，取值都在 $$[0,1]$$ 区间上；其中每条折线代表一朵花。从图中我们可以看出，setosa 这种花的花瓣较小，而另外两种花的花瓣都较大；但 setosa 的花萼特别宽；versicolor 这种花总体来说和 virginica 比较像，只是相对来说尺寸小一些。另外，从图中我们也可以看出相邻变量之间的正向或负向相关关系，例如，由于花萼宽（Sepal.Width）和花瓣长（Petal.Length）两个变量间的线段有很多相交，说明这两个变量有负相关关系，而对于 setosa 这类花来说，由于花萼长和花萼宽这两个变量之间的线段都朝向相同的方向，因此这两个变量在 setosa 这类花中有正向相关关系。这些结论都可以用数值方式验证：

# 第 2、3 列的相关系数（花萼宽和花瓣长）
cor(iris[, 2:3])
##              Sepal.Width Petal.Length
## Sepal.Width    1.0000000   -0.4284401
## Petal.Length  -0.4284401    1.0000000
# 第 1、2 列在 setosa 下的相关系数（先对数据取子集）
cor(subset(iris, Species == "setosa", 1:2))
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000   0.7425467
## Sepal.Width     0.7425467   1.0000000

## 4.37 调和曲线图

$$$f_{i}(t)=\begin{cases} \frac{X_{i,1}}{\sqrt{2}}+X_{i,2}\sin(t)+X_{i,3}\cos(t)+\cdots\\ \qquad+X_{i,p-1}\sin(\frac{p-1}{2}t)+X_{i,p}\cos(\frac{p-1}{2}t) & \text{若}p\text{为奇数}\\ \frac{X_{i,1}}{\sqrt{2}}+X_{i,2}\sin(t)+X_{i,3}\cos(t)+\cdots\\ \qquad+X_{i,p}\sin(\frac{p}{2}t) & \text{若}p\text{为偶数} \end{cases} \tag{4.17}$$$

1. 如果我们用 $$L_{2}$$ 范数来度量两条曲线之间的距离，那么得到的距离值正好是欧氏距离平方的 $$\pi$$ 倍，换句话说，两行观测之间的距离恰好可以表现为图中两条曲线之间的差距。这条性质使得我们可以直观地在图中观察聚类现象和离群点，因为聚类和离群点的概念都是基于距离的（距离的定义有多种，这里用欧氏距离的平方）。如果读者感兴趣，可以验证一下这个 $$L_{2}$$ 范数的结果： $\int_{-\pi}^{\pi}\left(f_{i}(t)-f_{j}(t)\right)^{2}dt=\pi\sum_{k=1}^{p}\left(X_{i,k}-X_{j,k}\right)^{2}$

2. 这个变换从一定程度上保持了线性性，即：若一个观测 $$X_{l}$$ 的所有数值都小于 $$X_{i}$$ 而大于 $$X_{j}$$，那么在调和曲线图上 $$X_{l}$$ 对应的曲线也位于 $$X_{i}$$$$X_{j}$$ 之间。这一点性质是非常明显的。

library(MSG)
usage(andrews_curve)
## andrews_curve(x, n = 101, type = "l", lty = 1, lwd = 1, pch = NA, xlab = "t",
##   ylab = "f(t)", ...)

iris.col <- vec2col(iris$Species) par(mfrow = c(2, 2)) par(mar = c(4, 4, 0.2, 0.2)) andrews_curve(iris[, 1:4], n = 50, col = iris.col, xlab = "$t$", ylab = "$f(t)$") andrews_curve(iris[, c(3, 4, 2, 1)], n = 50, col = iris.col, xlab = "$t$", ylab = "$f(t)$") andrews_curve(scale(iris[, 1:4]), n = 50, col = iris.col, xlab = "$t$", ylab = "$f(t)$") x <- andrews_curve(scale(trees), n = 50, xlab = "$t$", ylab = "$f(t)$") # 离群点是哪行数据？即哪行数据对应的 f(t) 会大于 4？ which(apply(x > 4, 1, any)) ## [1] 31 4.43 展示了两个多维数据的调和曲线图。左上图为原始鸢尾花数据的前四列，并根据鸢尾花的种类对曲线进行了颜色标注，可以看出，setosa 这一类花的曲线“拧”成了一股，说明这些花之间的距离比较近，也就是我们至少可以观察到这一类的聚类现象。右上图调整了四列的顺序，将花瓣长和花瓣宽两个变量放在了前面，这样做的原因我们后面再解释，这里我们首先可以看到，调整顺序后的鸢尾花调和曲线图发生了变化，曲线的振幅和频率都加大了，这跟数据各列的性质有关，鸢尾花数据有四列，我们可以写出曲线的函数表达式： $f_{i}(t)=\frac{X_{i,1}}{\sqrt{2}}+X_{i,2}\sin(t)+X_{i,3}\cos(t)+X_{i,4}\sin(2t)$ 我们知道，$$\sin(kt)$$ 中的 $$k$$ 越大，则三角函数曲线的频率越高，而每个三角函数前面的系数越大，则振幅越大。鸢尾花数据前四列的简单汇总信息如下： summary(iris[, 1:4]) ## 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 平均来说，第 1、2 列花萼比第 3、4 列花瓣的长宽数值要大，右上图中，由于大数值赋给了 $$\sin(2t)$$ 的系数，所以导致这一项相对占了主导地位，因此曲线频率加大；同时，由于曲线表达式的第一项理论上取值无界（后面所有项都受三角函数取值界限限制），所以它往往很大影响了曲线的振幅，右图中将小值给了第一项，所以振幅相对左上图加大了。右上图的另一个特征是，setosa 类的曲线几乎一致“拧”在所有曲线之下，这也使得观察聚类更方便。为什么要把第 3、4 列调到前面来呢？主要是因为鸢尾花的四个属性中，用花瓣长宽最易区分各类花，这一点，我们可以用分类树验证一下： library(rpart) print(rpart(Species ~ ., iris), digits = 2) ## n= 150 ## ## node), split, n, loss, yval, (yprob) ## * denotes terminal node ## ## 1) root 150 100 setosa (0.333 0.333 0.333) ## 2) Petal.Length< 2.5 50 0 setosa (1.000 0.000 0.000) * ## 3) Petal.Length>=2.5 100 50 versicolor (0.000 0.500 0.500) ## 6) Petal.Width< 1.8 54 5 versicolor (0.000 0.907 0.093) * ## 7) Petal.Width>=1.8 46 1 virginica (0.000 0.022 0.978) * 在数据各列的数量级都差不多的时候， (4.17) 式越靠前的变量越能控制曲线的位置，越靠后的变量越能控制曲线的波动，而我们通常并不关心调和曲线图中的波动，仅仅观察曲线是否“拧成股”，所以在安排变量顺序的时候，通常把对聚类有重要贡献的变量放在前面，这样调和曲线图对读者来说才能有更好的可读性。 4.43 左下图是将鸢尾花数据标准化之后的调和曲线图，由于数据的数量级更加靠近，所以曲线之间的差异也更明显，导致各股曲线更加松散。右下图是标准化之后的黑莓树数据调和曲线图，该数据名为 trees，在 datasets 包中，数据包含 3 列：周长、高和体积；从图中我们可以观察到有一条曲线明显与其它曲线“不合群”，因此可以初步判断它是一个离群点，由于这条曲线的数值能够超过 4，我们可以用函数返回值来查找它对应的行（事实是第 31 行），代码见图 4.43 ## 4.38 二维箱线图 4.3 小节中我们介绍了普通的箱线图，即用箱线表示一维数据的各个分位数，在二维情况下，我们可以用类似的思想画二维箱线图。二维箱线图又名袋图（Bag Plot），它由 提出。二维箱线图的做法是从数据的中心向外，逐渐用凸包多边形将散点图中的点包起来，直到包到一半的数据点，此时的凸包相当于普通箱线图中的箱子，然后再向外包到所有数据点。二维箱线图的基本构成就是一个中心和两个多边形，它们能粗略描述数据的二维分布情况。 R 中 aplpack 提供了一个函数 bagplot() 可以用来画二维箱线图，其用法如下： library(aplpack, warn.conflicts = FALSE) usage(bagplot) ## bagplot(x, y, factor = 3, na.rm = FALSE, approx.limit = 300, ## show.outlier = TRUE, show.whiskers = TRUE, show.looppoints = TRUE, ## show.bagpoints = TRUE, show.loophull = TRUE, show.baghull = TRUE, ## create.plot = TRUE, add = FALSE, pch = 16, cex = 0.4, dkmethod = 2, ## precision = 1, verbose = FALSE, debug.plots = "no", col.loophull = "#aaccff", ## col.looppoints = "#3355ff", col.baghull = "#7799ff", ## col.bagpoints = "#000088", transparency = FALSE, show.center = TRUE, ...) 其中 xy 分别是横纵坐标轴上的数据向量，也可以直接提供一个 2 列的矩阵或数据框； factor 类似 boxplot() 中的 range 参数，用来定义离群点，取值越大，则离群点越少（数据点离中心的距离可以越远）； approx.limit 界定了大数据的样本量，如果原始数据的样本量超过这个数字，则随机抽取 approx.limit 个数据点用作二维箱线图的计算； dkmethod 取值 1 或 2，决定用哪种方法计算袋子的范围，取值 2 计算更精确；其它设置颜色或形状的参数此处略去。 # 数据来自 rpart 包中的 car.test.frame par(mfrow = c(1, 2), mar = c(4, 4, 0.1, 0.1), bty = "l") car.wd <- with(car.test.frame, cbind(Weight, Disp.)) bagplot(car.wd, cex = 0.9) box() library(KernSmooth) # 二维核密度估计并画等高线图 est <- bkde2D(car.wd, apply(car.wd, 2, dpik)) with(est, { contour(x1, x2, fhat, xlab = "Weight", ylab = "Disp.", xlim = range(car.wd[, 1]), ylim = range(car.wd[, 2]) ) points(car.wd, pch = 16, cex = 0.9) }) 4.44 左图展示了 Consumer Reports 中汽车数据中车重 Weight 和气缸排量 Disp. 两个变量的二维箱线图，图中心有一个米字型的点，代表二维数据的中位数，内层深色袋子包含了一半的数据点（约 30 个），外层袋子本应包含所有数据点，但由于默认设置有 5 个离群点的存在，所以只包含了 55 个数据点，读者可以将 factor 参数设置为 6，便可包含所有数据点。作为对比，右图从二维核密度估计的角度画出了密度值的等高图，这两幅图可以看作是相通的，左图里层袋子就是右图中密度曲面5 从中心向外积分数值到 0.5 时对应的区域。 ## 4.39 延伸与小结 到目前为止，我们已经用了 38 小节的篇幅介绍了所有 graphics 包中的图形以及部分附加包中的图形，基于以下两个考虑，我们将在本小节简略介绍一下其它包中的其它图形之后结束本章： 1. R 附加包数目太多：截至 2019 年 8 月 25 日，CRAN 上的 R 包数量已经接近 14900 个，其中很多都包含作图的函数，如果我们这样继续列举下去，本书将永远没有尽头；有一个叫“R Graphical Manual”的网站，6基于 R 的附加包中的例子画出了所有图形，在 R 2.9.0 的 1877 个附加包的示例代码中，一共生成了 21924 幅图形，尽管有很多示例代码可以同时生成多幅图形，但两万多幅图形背后包含的图形种类很可能成百上千，这样的数目是本书无法承载的。 1. 很多 R 函数的用法都比较类似：一方面，有些附加包为了便于用户使用，只是扩展了泛型函数 plot()，这样用户只需要在建完特定的模型之后 plot() 相应的对象即可，例如 MASS 包中的岭回归；另一方面，对于基础图形而言，很多参数都是通用的，它们的意义通常比较固定，例如 col 参数通常表示主要图形元素的颜色等，这样，即使不看帮助文档，大致也能猜到用法，所以我们没有必要把每一个作图函数的用法都详细介绍一遍。 由于 R 拥有大量的附加包，图形种类也被极大扩展，R 的初级用户往往会提出这样一个自然的问题： 我怎么知道有哪些作图的包？或者说怎样找到合适我的包呢？ 这个问题并不容易回答，它需要结合具体问题来看。从一般原则来说，用户可以浏览一下 CRAN 关于图形的分类列表（Task View）：https://cran.r-project.org/web/views/Graphics.html，这是一个很好的导航页面，它总结了一些有用的包的功能。对于具体的问题，也可以广泛寻求网络上的帮助和指导（参见 2.2 小节）。 本小节选取了 plotrix 包中的两个作图函数作为代表，让读者了解基于 R 的基础图形系统的扩展可能性。 是在 R News 上关于这个包的一篇介绍性文章，文章标题将这个包描述为“（它）位于 R 的红灯区”，原因是它看起来触犯了一些作图的常见规则，比如它扩展了饼图（甚至还有 3D 饼图这种统计学家很不齿的图形！）、介绍了如何将坐标轴截断、可以画像 Excel 那样的渐变色条形图等等，这里本作者也不知该包的作者究竟是开玩笑还是认真的，总之这个包的确显示了 R 基础图形系统的扩展性，但一些扩展的方向需要我们三思，例如截断坐标轴的作法就是 所反对的。下面我们分别介绍风向图和浮动饼图。 library(plotrix) # 直接取自 oz.windrose() 函数的例子 windagg <- matrix(c( 8, 0, 0, 0, 0, 0, 0, 0, 4, 6, 2, 1, 6, 3, 0, 4, 2, 8, 5, 3, 5, 2, 1, 1, 5, 5, 2, 4, 1, 4, 1, 2, 1, 2, 4, 0, 3, 1, 3, 1 ), nrow = 5, byrow = TRUE) par(mar = c(0.5, 0.5, 4.0, 0.5)) oz.windrose(windagg) 风向图（Wind Rose）通常用来展示东南西北甚至更多方向上的风力强度和频数，它是观察风力在各个方向上分布的一种直观办法，以图 4.45 为例：它用八根指针表示东、东南、南、西南等八个方向上的风，在每一根指针上对各等风速进行频数汇总，最后处以观察时间内的刮风次数（包括无风次数），就得到了各种风速的频率，分别用指针的一小节来表示，也就是说，指针的某一小节越长，表示该风速水平上的刮风次数越多。中心的圆圈表示无风的频率，之所以要用圆圈，原因很简单，因为无风是没有风向的，无法用指针表达；无风的总频率被分配到了八个方向上，所以我们不能简单看圆圈的大小来决定无风频率，而是要将圆圈大小乘以方向的数目（这里为 8），当然这种计算对读者来说太麻烦，所以图的左边直接以文本标注的形式标出了无风的总频率。该风向图的模板取自澳大利亚气象局。 在一些与气象有关的部门，风向图有其特殊用处，例如机场需要参考当地的风向图来决定应该以什么样的方向修建飞机跑道，由于顺风对飞机飞行不利，所以飞机在机场起飞和降落最好是能逆风，这样能使飞机在较短距离升空和降落。若图 4.45 是一个真实的机场的长期风向数据，则机场跑道可能考虑东北-西南方向或者东西方向，前者的刮风频率较高，但后者的强度较大。作者对机场建设了解甚少，这里只是介绍一种考虑的可能性。另外，根据作者目前掌握的知识，当今的飞机似乎受风向影响已经大大减小了。 我们在风向图中看到一项数据处理的缺陷：风向数据本来应该是连续的，但在画图的时候无论是风向还是风速都被离散化了（也许风向在观测时就只能是离散的），这种离散化对读者理解来说有利，但它毕竟是损失了数据信息。若有气象工作者能看到这里，作者提议考虑将风速的图示表达方式换成小提琴图（4.33 小节），或者至少箱线图（4.3 小节），除非现在这种人工分组的方式有任何坚实的科学基础。 demo("float_pie", package = "MSG") 浮动饼图是 plotrix 包对饼图的一个扩展，它可以将饼图作为一个基本的图形元素安排在任意位置，即：一幅图中可以有多幅饼图分散在不同位置。尽管我们一直批评饼图在表达数据上的弱势，但在某些极端情况下，饼图也未必糟糕得无可救药，试想如果一幅饼图只分为两份或三份，也许我们可以准确识别各个角度的大小。另外，有时候我们也许要比较多组比例的分布情况，例如在空间中表示各个地点的某个变量构成（性别比例、年龄结构等），此时在多个位置分别画饼图也许能让我们很容易比较各个位置上的比例构成。图 4.46 展示了一个例子，改编自 floating.pie() 函数自身的例子，我们假设要检验一批产品的质量，若“黑色成分”超过 50 % 则表明质量不合格；图中我们将不合格的产品“沉”了下去，这种形象表达或许更能引起读图者的注意。 至此，除了 4.36 小节，我们所介绍的图形几乎都是基于 R 的基础图形系统构造的，即：用点、线、颜色、形状等基本元素来构造完整的图形。这种方式的确给了我们极大的自由，让我们可以随意控制图形，而且很多时候我们也不必用基本图形元素来自行构造图形，因为 R 及其附加包已经提供了太多相对完善的高层作图函数，这些函数应该可以涵盖相当一部分用户的需求。然而我们在最后需要提醒读者的是，在 R 的基础图形系统之外还存在好几种选择，我们的思维不必局限在如此原始的画图方式上。R 基础图形系统的缺点至少包括： 1. 图形元素一旦画出来就不可编辑，如果需要更改图形，只能重画整幅图形（实际上这个缺点并不严重，因为修改代码和重新运行代码对某些代码编辑器来说非常方便） 2. 作图功能不够自动，比如我们常常需要根据某个分类变量给相应的图形元素赋予颜色，这时候只能将分类变量人工转变为颜色向量，事后还要手工添加图例，再比如添加 LOWESS 光滑曲线或回归直线及其置信区间，用户需要经过繁琐的过程才能完成，而这些作图的常规任务应该能尽量自动完成，而不要让用户操心每一个细节； 3. 图形系统几乎没有交互功能，用户很难与图形设备交互，比如用鼠标选取一部分点并高亮之（附录 B.5 小节介绍了非常简单的交互功能）； 4. 有些细节设置不够合理，例如点的样式默认为空心圆圈（pch = 1），尽管空心圆圈在某些场合下有其特殊优点，但一般说来这种设置没有足够的视觉冲击力，点在图中不够明显，用 的话说，就是“数据不能突出来”； 接下来的一章，我们将介绍不同的图形系统来弥补 R 基础图形系统的不足。 ## 4.40 思考与练习 1. 在第 4.30 节中我们了解了如何画 QQ 图，与之对应的还有一种图形叫 PP 图（Probability-Probability），它也是一种检验数据分布是否和理论分布吻合的图形工具，原理和 QQ 图类似：对数据的实际概率分布值和理论概率分布值作散点图即可（也可以选择性地添加一条直线）。编写函数画出 geyser$waiting 的 PP 图（理论分布选择正态分布，均值和标准差用矩估计获得），并评价该数据的正态性。本题源于 COS 论坛帖子：https://d.cosx.org/d/18521

2. 用 QQ 图评价数据正态性应该注意什么问题？模拟一些正态分布的数据，画出 QQ 图，看看如果数据真的来自正态分布，QQ 图看起来是什么样的。或者参考 animation 包中的 sim.qqnorm() 函数。

3. 如果箱线图中出现大量的占据范围很大的离群点，则会导致箱子的主体部分被压得很扁，此时我们很难看清分位数的位置。请问对于这种情况有什么好的解决办法？

4. 条形图的横坐标通常是无序的，人们往往根据数据的原始顺序画条形图；比如画各省市的人口总数时，我们总是看到北京在第一位，这可能是因为统计局数据总是把北京放在第一行。条形图横坐标的顺序选取对读图是否有影响？或者说我们应该怎样安排条形图的横坐标？提示：参考 6.2.2 小节。

5. 尽管一元函数曲线图和展示数据似乎没什么关系，但它在优化一元目标函数或者对一元函数求根时通常有帮助。例如我们用 uniroot() 求根时需要提供根的大致区间，请参考帮助文档说明一元函数曲线图在这里有什么用处。

6. 基于第 4.15 小节介绍的矩阵图，编写一个画平行坐标图的函数，参数主要包括一个数据框（可以选择性地包括其它修饰性参数，例如数据标准化的方法等），函数的主体部分只有一行代码，形式如下：

parcoords <- function(x, ...) {
# 如何处理 x？
matplot(..., type = "l", lty = 1, pch = NA)
}
# 测试代码
parcoords(iris)

提示：可以考虑 col() 函数，可能需要转置 t() ；想清楚我们是用怎样的两个矩阵去画线。

7. 在信息可视化（Information Visualization）领域，树图和标签云都非常流行，实际上这些图形从统计学角度来说表达的信息非常简单：它们表达的只是数字大小。树图为 Treemap，不是 4.32 小节提到的分类与回归树，但思想类似，也是递归分割，如图 4.47 是本作者操作系统中的 R 包文件大小树图。树图主要是用矩形大小代表数字大小，纵横交替分割一个大矩形为小单元。图 4.47 中最大的灰色矩形代表所有 R 包的大小，在 ~/R/x86_64-pc-linux-gnu-library/2.12/ 文件夹下有若干 R 包，每个包都有自己的大小，其中最大的是 mapdata，其次是 RGtk2，然后是 Rcpp，等等。纵向划分的长矩形表示一个个包（面积和包的大小成比例），然后在这些矩形内部再横向划分小矩形分别表示子文件夹的大小，若子文件夹下还有子文件夹，那么继续纵横划分。树图同时表达了嵌套关系和数值大小，看起来一目了然，通常最大的矩形能最先吸引人的注意力。例如，若我们想清理磁盘上的文件，那么就会考虑究竟是哪些文件占用了很大的空间，此时树图就是很好的可视化方法（本作者正是用这样的办法清理自己的磁盘空间的）。图 4.47 显示的是 Ubuntu 系统下的 Disk Usage Analyzer，Windows 下也有类似的软件如 WinDirStat。

标签云就更简单了：将一些文本标签按一定顺序排在平面上，文本大小和某个数值成比例，这样最大的文字就能最先吸引我们的注意力。例如我们计算一篇文章中单词出现的频数，用频数大小来决定单词的大小，这样高频出现的词在图中一眼就能看出来，如图 4.48 是作者英文博客（https://yihui.org/en/）利用 Wordle 生成的标签云。显然，animation 和 Sweave 两个词字号最大，这说明这两个词在博客中出现频率最高，这不难解释：因为作者一直开发 animation 包，并研究着 Sweave。

本书大部分内容都是关于数据可视化的，统计学的成分比较大。请结合这里给出的两幅图形思考，信息可视化和数据可视化的区别是什么？它们各有什么优势以及如何结合它们的优势？

### 参考文献

Adler, Daniel, and S. Thomas Kelly. 2021. Vioplot: Violin Plot. https://github.com/TomKellyGenetics/vioplot.
Adler, Daniel, and Duncan Murdoch. 2021. Rgl: 3d Visualization Using OpenGL. https://CRAN.R-project.org/package=rgl.
Andrews, D. F. 1972. “Plots of High Dimensional Data.” Biometrics 28 (1): 125–36. https://doi.org/10.2307/2528964.
Breiman, Friedman, Olshen, and Stone. 1984. Classification and Regression Trees. Wadsworth.
Brownrigg, Ray. 2018. Maps: Draw Geographical Maps. https://CRAN.R-project.org/package=maps.
Chambers, J. M., W. S. Cleveland, B. Kleiner, and P. A. Tukey. 1983. Graphical Methods for Data Analysis. Wadsworth; Brooks/Cole.
Chernoff, Herman. 1973. “Using Faces to Represent Points in k-Dimensional Space Graphically.” Journal of the American Statistical Association 68 (342): 361–68. https://doi.org/10.2307/2284077.
Cleveland, William S. 1979. “Robust Locally Weighted Regression and Smoothing Scatterplots.” Journal of American Statistical Association 74 (368): 829–36. https://doi.org/10.2307/2286407.
———. 1985. The Elements of Graphing Data. Monterey, CA: Wadsworth.
———. 1993. Visualizing Data. Hobart Press.
Cohen, Ayala. 1980. “On the Graphical Display of the Significant Components in Two-Way Contingency Tables.” Communications in Statistics - Theory and Methods 9 (10): 1025–41. https://doi.org/10.1080/03610928008827940.
Friendly, M. 1992. “Graphical Methods for Categorical Data.” In SAS User Group International Conference Proceedings, 17:190–200.
Friendly, Michael. 1994. “A Fourfold Display for 2 by 2 by $$k$$ Tables.” 217. York University, Psychology Department. http://www.math.yorku.ca/SCS/Papers/4fold/4fold.ps.gz.
Hintze, J. L., and R. D. Nelson. 1998. “Violin Plots: A Box Plot-Density Trace Synergism.” The American Statistician 52 (2): 181–84. https://doi.org/10.1080/00031305.1998.10480559.
Hofmann, H., and M Theus. 2005. Interactive Graphics for Visualizing Conditional Distributions.
Hothorn, Torsten, and Brian S. Everitt. 2014. A Handbook of Statistical Analyses Using r. 3rd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://CRAN.R-project.org/package=HSAUR3.
Inselberg, Alfred. 2009. Parallel Coordinates: Visual Multidimensional Geometry and Its Applications. Springer-Verlag New York. https://doi.org/10.1007/978-0-387-68628-8.
Kaplan, E. L., and P. Meier. 1958. “Nonparametric Estimation from Incomplete Observations.” Journal of the American Statistical Association 53 (282): 457–81. https://doi.org/10.2307/2281868.
Lemon, Jim. 2006. plotrix: A Package in the Red Light District of r.” R News 6 (4): 8–12. https://cran.r-project.org/doc/Rnews/Rnews_2006-4.pdf.
Ligges, Uwe, and Martin Maechler. 2003. “Scatterplot3d - an r Package for Visualizing Multivariate Data.” Journal of Statistical Software 8 (11): 1–20. https://www.jstatsoft.org/v08/i11.
Meyer, David, Achim Zeileis, and Kurt Hornik. 2006. “The Strucplot Framework: Visualizing Multi-Way Contingency Tables with Vcd.” Journal of Statistical Software 17 (3): 1–48. http://www.jstatsoft.org/v17/i03.
———. 2020. Vcd: Visualizing Categorical Data. https://CRAN.R-project.org/package=vcd.
Robert McGill, John W. Tukey, and Wayne A. Larsen. 1978. “Variations of Box Plots.” The American Statistician 32 (1): 12–16. https://doi.org/10.2307/2683468.
Rousseeuw, P. J., I. Ruts, and J. W. Tukey. 1999. “The Bagplot: A Bivariate Boxplot.” The American Statistician 53 (4): 382–87. https://doi.org/10.1080/00031305.1999.10474494.
Sarkar, Deepayan. 2008. Lattice: Multivariate Data Visualization with r. New York, NY: Springer. http://lmdvr.r-forge.r-project.org.
Schloerke, Barret, Di Cook, Joseph Larmarange, Francois Briatte, Moritz Marbach, Edwin Thoen, Amos Elberg, and Jason Crowley. 2021. GGally: Extension to Ggplot2. https://CRAN.R-project.org/package=GGally.
Scott, David W. 2015. Multivariate Density Estimation: Theory, Practice and Visualization. 2nd ed. New York, NY: John Wiley & Sons, Inc. https://doi.org/10.1002/9780470316849.
Snow, Greg. 2020. TeachingDemos: Demonstrations for Teaching and Learning. https://CRAN.R-project.org/package=TeachingDemos.
Therneau, T. M., and P. M. Grambsch. 2000. Modeling Survival Data: Extending the Cox Model. New York, USA: Springer.
Therneau, Terry M. 2021. Survival: Survival Analysis. https://github.com/therneau/survival.
Therneau, Terry, and Beth Atkinson. 2019. Rpart: Recursive Partitioning and Regression Trees. https://CRAN.R-project.org/package=rpart.
Urbanek, Simon, and Tobias Wichtrey. 2018. Iplots: iPlots - Interactive Graphics for r. https://CRAN.R-project.org/package=iplots.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with s. 4th ed. New York, NY: Springer. https://www.stats.ox.ac.uk/pub/MASS4.
Wand, Matt. 2021. KernSmooth: Functions for Kernel Smoothing Supporting Wand & Jones (1995). https://CRAN.R-project.org/package=KernSmooth.
Wolf, Hans Peter. 2019. Aplpack: Another Plot Package: Bagplots, Iconplots, Summaryplots, Slider Functions and Others. http://www.wiwi.uni-bielefeld.de/lehrbereiche/statoekoinf/comet/wolf/wolf_aplpack.
Xie, Yihui. 2007. Visualization of Data and Statistical Models Using r. 第 3 届计量模型与中国·澳大利亚自由贸易区发展国际研讨会论文集. http://1t.click/ajpa.
———. 2020. MSG: Data and Functions for the Book Modern Statistical Graphics. https://github.com/yihui/MSG.
Zeileis, Achim, David Meyer, and Kurt Hornik. 2007. “Residual-Based Shadings for Visualizing (Conditional) Independence.” Journal of Computational and Graphical Statistics 16 (3): 507–25. https://www.jstor.org/stable/27594258.

1. 本书作者在爱荷华州立大学统计系读博士期间，每周统计图形小组有一次讨论，这则消息来自于作者的一位导师 Heike Hofmann 教授↩︎

2. 对于连续的自变量则称为 Spine Plot，离散自变量则为 Spinogram，这个单词是仿照 Histogram 而造。↩︎

3. 若因变量为分类变量，则预测值按照多数投票表决（majority vote）原则计算；若为数值变量，则按照叶节点上的样本均值预测。↩︎

4. 不要忘记：等高图实际上是三维图形！请想象这个“曲面”的形状。↩︎