RColorBrewer
文档按数据类型选预设色盘。
章节概述
数据可视化帮助观察分布、差异、变量关系与时间变化,也让读者能够核对研究论述。本章从为什么看图、怎样读图出发,再用
ggplot2 的数据与映射 → 绘制图层 →
标度与主题建立构图框架。
九类常用静态图依次说明适用问题、数据准备、最小代码、读图与进阶表达;动态与交互示例展示年份动画、切换、悬停和三维视角。学术制图还需检查选图、尺度、单位、坐标、图例和标注,并通过配色、字体与布局调整,提高图形的清晰度与一致性。
基本概念
数据可视化 (data visualisation) 把变量的取值对应到图上的位置、长度、颜色或形状。例如,散点图可以用一个点表示一条记录,用横纵坐标表示两个数值变量,用颜色区分组别。这种对应关系称为视觉映射。
统计表便于查找具体数值,汇总指标便于概括样本;图形便于观察分布、组别差异和变量关系。读图时仍须确认记录代表谁、坐标轴使用什么单位,以及图中出现的模式是否受个别观测影响。
看图主要为了做三件事:检查数据形态、比较差异与变化、说明研究发现。
检查数据形态 分布图可以显示数值集中在哪些范围;散点图可以显示两个变量的关系是否接近直线、是否弯曲,以及有无明显偏离多数点的观测。图上的异常是进一步核对原始记录的线索。缺失值通常不会自动出现在常规散点图中,需要另行统计或绘图检查。
统计学家 Frank Anscombe 在 1973 年给出四组各 11 条观测。四组的 \(x\) 均值都是 9,\(y\) 均值约为 7.50;\(x\) 的样本方差都是 11,\(y\) 的样本方差约为 4.12—4.13,相关系数约为 0.816,线性拟合结果都接近 \(\hat y = 3 + 0.5x\)。这些统计摘要十分接近,散点形态却明显不同。
这里的方差概括数值有多分散,相关系数概括线性关系的方向与强弱;\(\hat y\) 是拟合线对 \(y\) 的预测值。暂不需要手算这些指标,先看四幅散点图。
安斯库姆四重奏:相近统计摘要下的四种散点形态;蓝点为原始观测,橘红线为各组线性拟合
两个读图术语。 离群点 (outlier) 是明显偏离多数观测的点。高杠杆点 (high-leverage point) 的横坐标远离多数观测,可能明显影响回归线。发现这类点时,应先核对记录和分析口径,不能仅凭图形将其删除。
四组相近的摘要对应了不同的数据形态。散点图提示了曲线关系和特殊观测;线性模型是否合适、个别观测影响有多大,还需结合后续检查判断。
比较差异与变化 并列图形让组别差异更容易比较,按时间排列的点和线让变化方向更清楚;误差棒与阴影带还可以提示估计存在不确定性。下图示意这三种用途。
比较组别、观察趋势与表达不确定性的概念示意
说明研究发现 图表可以让组别差异和变化方向更容易查找,也便于读者对照正文中的论述。图形的说服力仍取决于数据质量、分析方法以及标题、单位、图例和来源是否交代清楚。
图形在研究结果沟通中的用途示意
第二章使用过 gapminder
数据:一行记录一个国家或地区在某一年的情况。本章接着提出一个具体问题:只看
2007 年的记录,各国或地区的人均 GDP
与出生时预期寿命呈现什么关系? 该教学数据覆盖 142
个国家或地区,不能视为完整的全球人口数据。
散点图可以让一个点对应一个国家或地区;横轴放人均
GDP,纵轴放出生时预期寿命,再用颜色区分大洲。预期寿命以岁为单位。随后观察点群的总体方向、大洲之间的分布,以及远离多数点的观测。2.2
节将用 ggplot2 绘制这张图。
读图时可依次核对三个问题:
即使点群呈现正向关系,也只能说明这份样本中的两项指标有关联。要判断人均 GDP 是否影响预期寿命,还需考虑其他因素和研究设计。
出图规范将在 4.0 专门讨论。 例如,用柱长或条长表示非负数量时,数值轴通常从 0 开始;散点图可按数据范围设置坐标,但需清楚标出刻度和变换。这里先把注意力放在数据口径与读图判断上。
工具基础
ggplot2 借鉴了统计学家 Leland Wilkinson
在《图形语法》(The Grammar of
Graphics)中提出的思路。一张图可以从基本设置开始,按需要逐步组合:
下图用 R 自带的汽车数据 mtcars
演示这种组合方式:wt 是车重(千磅),mpg
是每加仑汽油可行驶的英里数。先看图中的完整代码,2.1
节再按顺序拆开阅读。
ggplot2 构图示意:数据与变量映射建立绘图基础,点和拟合线形成图层,主题调整外观
上方示意图的底部给出了一段完整代码。沿着箭头阅读,可以看到每一部分怎样改变前一步的图形:
mtcars %>% ggplot(aes(x = wt, y = mpg)) 将
mtcars 交给 ggplot();aes() 把
wt 指定为横轴、mpg
指定为纵轴。数据中一行对应一款车型的记录。到这一步,图形有了数据和位置规则,还没有画出点。+ geom_point()
按照刚才的横纵轴规则,把每款车型画成一个点。只运行到这一步,已经能看到车重与燃油效率的分布。+ stat_smooth(method = "lm")
在点上叠加一条线性拟合线,用直线概括两项指标在这份数据中的整体方向;线周围的浅色带提示拟合结果仍有不确定性。此处只需识别它在图中的作用。+ theme_minimal()
将背景和网格整理为简洁样式,不改变原始数据,也不改变横纵轴代表的变量。前两步已经能画出基础散点图;拟合线与主题设置可按表达需要继续添加。2.2 节将用同样的组合方式,换一份社会科学教学数据实际绘图。
ggplot2 与
gapminder
基础 · 代码实战
本节用 gapminder 的 2007 年记录回答 1.2
节提出的问题:各国或地区的人均 GDP
与出生时预期寿命呈现怎样的关系?
每个点代表一个国家或地区,颜色区分数据中的大洲分组。先用少量代码把数据、变量映射和点图层接起来。
数据口径:Gapminder 教学数据。
本章使用 gapminder::gapminder 的历史教学子集:包含 142
个国家或地区,1952—2007 年每隔五年记录一次,共 12
个时点。一行是一个国家或地区在某一年的记录。该数据并非持续更新的完整全球数据。
| 变量 | 含义与单位 |
|---|---|
country、continent |
国家/地区及其在数据中的大洲分组 |
year |
观测年份,1952、1957、……、2007 |
lifeExp |
出生时预期寿命,单位为岁;不是居民平均年龄或平均死亡年龄 |
pop |
人口数量,单位为人 |
gdpPercap |
人均 GDP;数据集文档注明以经通胀调整的美元计 |
后文提到的“各国”“各洲”,均指这份教学数据所覆盖的记录。
运行前检查。 本节使用 tidyverse 与
gapminder。首次使用时,分别在 Console 运行
install.packages("tidyverse") 和
install.packages("gapminder");安装完成后,在当前 R
会话中加载:
library(tidyverse)
library(gapminder)
data_2007 <- gapminder %>%
filter(year == 2007)
ggplot(data_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) +
geom_point() +
labs(x = "人均 GDP", y = "出生时预期寿命(岁)", color = "大洲")逻辑对应
ggplot(data_2007, ...) 指定绘图记录;aes()
将人均
GDP、预期寿命和大洲分别交给横轴、纵轴和颜色;geom_point()
把一条记录画成一个点。labs()
只补充读图所需的变量名称和单位,不改变数据或图层。
回看 1.2 的读图问题。 图中每个点是教学样本中一个国家或地区在 2007 年的记录。多数点集中在横轴较低的区间,人均 GDP 较高的点大多位于预期寿命较高的区域;各洲的分布也有重叠。图形呈现的是这份样本中的同年关联,不能据此断言人均 GDP 的变化导致预期寿命变化。
ggplot2
实战进阶
进阶 · 出图展示
下面沿用 2007 年的 142 条记录,把人口规模加入气泡大小,并叠加一条基于全部样本的平滑趋势。图中还用简短标签标出四个国家,便于定位。先看成图,再按需展开完整代码;现阶段无需逐行掌握进阶设置。
为什么改用对数横轴? 在 2.2
的线性横轴上,多数点挤在较低数值区间。对数刻度让 500 到 5,000 与 5,000
到 50,000
这两段相同倍数的变化占据相同距离,便于观察低值区间的差异;刻度标签仍显示原始的人均
GDP。这里所有人均 GDP 都大于
0,符合对数变换的要求。scale_x_log10() 也会使下方的 LOESS
曲线基于 log10(人均 GDP) 拟合。
2007 年教学样本的人均 GDP、出生时预期寿命、人口规模和大洲分组
读图核对。 这些记录总体上沿人均 GDP 较高、预期寿命也较高的方向分布,但同一 GDP 区间仍有明显差异。趋势线由全部 142 条记录共同估计,不代表每个大洲各自的走势;图形也不能单独证明因果关系。
进阶代码提示: 气泡大小、趋势线、标签、配色和排版增加了设置,但构图仍沿用“数据与映射 → 绘制图层 → 标度 → 主题”的基本框架。
完整代码可复制运行。运行前需安装
ggrepel(调整标签位置)与 ragg(Knit
时绘图);首次使用时,在 Console 运行
install.packages(c("ggrepel", "ragg"))。图中只调用本机已有字体,不下载网络字体。绘图对象
plot_advanced 会留在当前 R
会话中,调整颜色或标签后可重新运行。
# 0. 环境与本地字体
library(tidyverse)
library(gapminder)
library(ggrepel)
font_priority <- c(
"Microsoft YaHei", "Microsoft YaHei UI", # Windows
"PingFang SC", "Hiragino Sans GB", "Heiti SC", # macOS
"Noto Sans CJK SC", "Source Han Sans SC" # 已安装的其他中文字体
)
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
installed_fonts <- unique(systemfonts::system_fonts()$family)
font_candidates <- font_priority[font_priority %in% installed_fonts]
supports_chinese <- function(family) {
tryCatch(
all(systemfonts::glyph_info(c("中", "文", "R"), family = family)$index > 0),
error = function(e) FALSE
)
}
font_candidates <- font_candidates[
vapply(font_candidates, supports_chinese, logical(1))
]
if (length(font_candidates)) font_family <- font_candidates[[1]]
}
# 1. 数据准备:大气泡先画、小气泡后画;四个名称仅作定位标签
data_2007 <- gapminder %>%
filter(year == 2007) %>%
arrange(desc(pop))
stopifnot(nrow(data_2007) == 142L, all(data_2007$gdpPercap > 0))
label_data <- data_2007 %>%
filter(as.character(country) %in% c("China", "India", "United States", "Japan")) %>%
mutate(country_label = recode(as.character(country),
China = "中国", India = "印度", `United States` = "美国", Japan = "日本"
))
continent_cols <- c(
Africa = "#B86E58", Americas = "#4D79A8", Asia = "#288E86",
Europe = "#8173A8", Oceania = "#B28E45"
)
trend_label <- "全样本 LOESS"
# 2. 绘图:先画趋势,再画气泡;之后调整标度、标签、图例与主题
plot_advanced <- ggplot(
data_2007, aes(x = gdpPercap, y = lifeExp)
) +
# 图层 A:全样本平滑趋势与拟合均值的近似 95% 置信带
geom_smooth(
aes(linetype = trend_label),
method = "loess", formula = y ~ x, se = TRUE,
colour = "#43536B", fill = "#A9BBCB",
alpha = 0.19, linewidth = 0.95,
show.legend = NA
) +
# 图层 B:填充色映射大洲,显示大小映射人口(单位:亿人)
geom_point(
aes(size = pop / 1e8, fill = continent),
shape = 21, colour = "#3B4C5D", stroke = 0.18,
alpha = 0.72
) +
# 标度:log10 横轴压缩高值区间;具名色值固定大洲配色
scale_x_log10(
breaks = c(300, 1000, 3000, 10000, 30000),
labels = scales::label_number(big.mark = ","),
minor_breaks = NULL,
expand = expansion(mult = c(0.025, 0.05))
) +
scale_y_continuous(
breaks = seq(40, 80, by = 10),
expand = expansion(mult = c(0.055, 0.06))
) +
scale_fill_manual(
values = continent_cols, breaks = names(continent_cols),
labels = c("非洲", "美洲", "亚洲", "欧洲", "大洋洲"), name = "大洲"
) +
scale_size(
name = "人口(亿人)", range = c(1.7, 11.5),
limits = c(0, 14), breaks = c(0.5, 2, 5, 10)
) +
scale_linetype_manual(
name = "总体趋势", values = setNames("dashed", trend_label)
) +
# 四个标签只是定位辅助;全部 142 个气泡仍在图中
ggrepel::geom_label_repel(
data = label_data, aes(label = country_label),
seed = 20260926, size = 3.15, family = font_family,
colour = "#263442", fill = scales::alpha("white", 0.93),
label.size = 0, label.padding = grid::unit(0.12, "lines"),
box.padding = 0.34, point.padding = 0.28, max.overlaps = Inf,
min.segment.length = 0, segment.size = 0.28,
segment.colour = "#8B99A6", show.legend = FALSE
) +
guides(
fill = guide_legend(
order = 1, override.aes = list(shape = 21, size = 4.2,
alpha = 1, colour = "#3B4C5D")
),
size = guide_legend(
order = 2, override.aes = list(shape = 21, fill = "#B8C5D0",
colour = "#43536B", alpha = 0.9)
),
linetype = guide_legend(
order = 3,
# 图中保持虚线;图例线段较短,改用实线便于辨认
override.aes = list(linetype = "solid", colour = "#43536B",
fill = "#A9BBCB", linewidth = 1.1)
)
) +
labs(
title = "人均 GDP 与预期寿命",
subtitle = "2007 年 · 142 个国家或地区 · 气泡大小表示人口规模",
x = "人均 GDP(对数尺度)",
y = "出生时预期寿命(岁)",
caption = "数据:gapminder。虚线为全样本 LOESS;浅色带为拟合均值的近似 95% 置信区间。"
) +
theme_minimal(base_size = 11.3, base_family = font_family) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_line(colour = "#EEF1F4", linewidth = 0.28),
panel.grid.major.y = element_line(colour = "#E4E9EE", linewidth = 0.34),
axis.line.x = element_line(colour = "#647180", linewidth = 0.38),
axis.line.y = element_line(colour = "#647180", linewidth = 0.38),
axis.ticks = element_line(colour = "#647180", linewidth = 0.32),
axis.ticks.length = grid::unit(2.4, "mm"),
axis.title = element_text(colour = "#314051", size = 10.5),
axis.text = element_text(colour = "#566575", size = 9.4),
plot.title = element_text(colour = "#243446", size = 18, face = "bold",
margin = margin(b = 4)),
plot.subtitle = element_text(colour = "#627181", size = 10.2,
margin = margin(b = 15)),
plot.caption = element_text(colour = "#718091", size = 8.3,
hjust = 0, margin = margin(t = 13)),
plot.title.position = "plot",
plot.caption.position = "plot",
legend.position = "right",
legend.box = "vertical",
legend.title = element_text(colour = "#304051", size = 9.7, face = "bold"),
legend.text = element_text(colour = "#445464", size = 8.8),
legend.key = element_rect(fill = "white", colour = NA),
legend.key.width = grid::unit(13, "mm"),
legend.spacing.y = grid::unit(3, "mm"),
plot.margin = margin(14, 18, 12, 14)
)
plot_advanced
深度解析
2.3 节的进阶示例延续基础散点图,在同一份 2007 年数据上增加人口大小、总体趋势与排版设置。此处先辨认各部分的职责;完整代码可在折叠区对照查看。
filter(year == 2007)
保留同一年的记录。ggplot(data_2007, aes(x = gdpPercap, y = lifeExp))
将横轴和纵轴设为所有图层共用的映射,因此趋势线和气泡都使用相同的两个变量。
geom_smooth() 先画出全部 142 条记录的一条 LOESS
平滑曲线,再由 geom_point() 画出各国或地区的气泡。LOESS
根据邻近点拟合平滑走势,无需预先指定一条直线。曲线周围的浅色带表示这条趋势线估计的不确定性(近似
95%
置信区间),并非单个国家的取值或未来结果的范围。气泡画在后面,因而显示在曲线之上。
气泡的 fill = continent 和 size = pop / 1e8
写在 aes()
中,分别让大洲决定填充色、让人口(亿人)决定气泡大小;边框颜色等固定样式写在
aes()
外。代码为五个大洲分别指定颜色,重绘时不会因记录顺序变化而换色。气泡略带透明度,大气泡先画、小气泡后画,使重叠处的记录更容易辨认。scale_x_log10()
改变横轴刻度,scale_size()
控制气泡显示尺寸;这些是标度设置,不是新的图层。
labs() 写明标题、单位和来源,guides()
调整图例符号,theme()
统一字体、留白和网格线。图中的总体趋势使用虚线;图例中的短线段改为实线,是为了在有限长度内更清楚地辨认这一图例项,两者指向同一条曲线。少数国家名称仅用于帮助定位,所有样本记录仍保留在图中。
区分三个位置。
ggplot()
指定数据与各图层共用的变量映射。geom_*()
增加真正画在图上的元素;某层独有的映射可写在该层的 aes()
中。scale_*()、labs() 与
theme()
分别调整变量如何显示、图中文字与整体外观。通常,随记录改变的属性写进
aes(),固定样式写在外面;本图为生成总体趋势图例,额外在
aes() 内映射了一个统一标签。3.0 节将进一步讨论散点图、折线图等图形分别适合回答什么问题。
本节依次认识九类常用图形。每一类都从实际问题出发,先准备数据,再用尽量少的绘图代码得到基础图;图出现后,再练习读图、观察哪些地方值得改进,并查看进阶成图。3.10—3.12 另设交互图、动画和 AI 辅助作图三个选学进阶模块,可按需要展开。
常用统计图形
本节选取九类常见图形,展示它们适用的数据、基础画法和进阶表达。多数示例使用
ggplot2 及其扩展包;雷达图的基础图使用
fmsb,进阶图使用 ggplot2。
点击卡片进入对应图形。
💡 学习建议:先掌握数据、映射与绘制内容如何组合,再留意各类图对数据形状和统计口径的不同要求。雷达图的两种画法会在该节单独说明。
交互图适合让读者主动改变观察角度。 下方四组八幅图两两并列:同一主题分别呈现原始观测、分布摘要或另一种视角。拖动、切换、悬停和旋转均可在本页操作;阅读时仍须核对坐标、色标与统计口径。
选用哪种工具? plotly
可把悬停、切换和动画嵌入本章的静态 HTML;shiny
可以让筛选条件触发 R 重新计算,需运行本地 app.R
或部署应用。gganimate 适合输出
GIF、视频等自动播放的动画,不能在图中悬停查询或筛选。若需将 Shiny
单独部署到静态站点,还可使用 shinylive
导出应用资源;本节不依赖它。
运行准备。 本页的八幅 Plotly 图已嵌入
HTML。复制代码到 RStudio 前,首次在 Console 运行
install.packages(c("plotly", "gapminder", "dplyr", "shiny", "bslib", "viridisLite"));绘图代码中的
library() 仍需随脚本运行。三维图使用 R 自带的
volcano 数据。
两图都使用 Gapminder 的五年间隔记录。左图逐年展示国家点云与当年全体样本的线性回归;右图把四个大洲样本的人口加权寿命连成逐步延长的轨迹。拖动滑块可分别查看横截面关系和长期变化。
A1 · 散点、回归线与 95% 置信带
点色区分大洲;悬停可读国家、人口及坐标。方程每帧重新拟合,浅色带是拟合均值的逐点 95% 置信区间,不能解释为单个国家的预测区间或因果效应。
library(dplyr)
library(gapminder)
library(plotly)
continent_cols <- c(
Africa = "#C35A68", Americas = "#3477A4", Asia = "#2F8D73",
Europe = "#8868A8", Oceania = "#B28B32"
)
continent_zh <- c(
Africa = "非洲", Americas = "美洲", Asia = "亚洲",
Europe = "欧洲", Oceania = "大洋洲"
)
# 每个国家每个年份保留一条记录;各洲颜色在全部帧中固定。
scatter_data <- gapminder %>%
mutate(
year = factor(year, levels = sort(unique(year))),
continent = factor(as.character(continent), levels = names(continent_cols)),
hover_label = sprintf(
"%s<br>%s 年 · %s<br>人均 GDP:%s 美元<br>预期寿命:%.1f 岁<br>人口:%.1f 百万人",
country, year, continent_zh[as.character(continent)],
format(round(gdpPercap), big.mark = ",", scientific = FALSE),
lifeExp, pop / 1e6
)
) %>%
arrange(continent, country, year)
# 每年单独拟合;置信带表示拟合均值的逐点 95% 置信区间。
fit_data <- scatter_data %>%
group_by(year) %>%
group_modify(~ {
fit <- lm(lifeExp ~ log10(gdpPercap), data = .x)
x_grid <- exp(seq(log(min(.x$gdpPercap)),
log(max(.x$gdpPercap)), length.out = 90))
ci <- predict(fit, newdata = data.frame(gdpPercap = x_grid),
interval = "confidence", level = 0.95)
tibble(gdpPercap = x_grid, fitted = unname(ci[, "fit"]),
lower = unname(ci[, "lwr"]), upper = unname(ci[, "upr"]))
}) %>%
ungroup()
fit_labels <- scatter_data %>%
group_by(year) %>%
group_modify(~ {
fit <- lm(lifeExp ~ log10(gdpPercap), data = .x)
coef <- unname(stats::coef(fit))
tibble(
gdpPercap = 60000, lifeExp = 26,
equation = sprintf("%s 年:ŷ = %.1f + %.1f × log₁₀(x)",
as.character(.y$year), coef[1], coef[2])
)
}) %>%
ungroup()
scatter_plot <- plotly::plot_ly(height = 570)
for (group_name in names(continent_cols)) {
group_data <- scatter_data %>% filter(continent == group_name)
scatter_plot <- scatter_plot %>%
plotly::add_markers(
data = group_data, x = ~gdpPercap, y = ~lifeExp,
frame = ~year, ids = ~country,
marker = list(color = unname(continent_cols[group_name]), size = 8.5,
opacity = 0.74,
line = list(color = "rgba(255,255,255,0.65)", width = 0.4)),
text = ~hover_label, hoverinfo = "text",
name = unname(continent_zh[group_name]), showlegend = TRUE,
inherit = FALSE
)
}
scatter_plot %>%
plotly::add_ribbons(
data = fit_data, x = ~gdpPercap, ymin = ~lower, ymax = ~upper,
frame = ~year, inherit = FALSE,
fillcolor = "rgba(93,150,181,0.20)",
line = list(color = "rgba(93,150,181,0)"),
name = "拟合均值的逐点 95% 置信区间",
hoverinfo = "skip", showlegend = TRUE
) %>%
plotly::add_lines(
data = fit_data, x = ~gdpPercap, y = ~fitted, frame = ~year,
inherit = FALSE, line = list(color = "#557F99", width = 3),
name = "当年样本回归线", hoverinfo = "skip", showlegend = TRUE
) %>%
plotly::add_text(
data = fit_labels, x = ~gdpPercap, y = ~lifeExp,
text = ~equation, frame = ~year, inherit = FALSE,
textposition = "top left", textfont = list(size = 12),
hoverinfo = "skip", showlegend = FALSE
) %>%
plotly::layout(
xaxis = list(title = "人均 GDP(美元,对数刻度)", type = "log", range = c(2.35, 5.1)),
yaxis = list(title = "出生时预期寿命(岁)", range = c(20, 88)),
legend = list(orientation = "v", x = 1.01, y = 1),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif"),
paper_bgcolor = "rgba(0,0,0,0)", plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 65, r = 185, b = 60, t = 15)
) %>%
plotly::animation_opts(900, easing = "linear", redraw = TRUE) %>%
plotly::animation_slider(currentvalue = list(prefix = "年份:")) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)A2 · 四洲人口加权寿命趋势
每个点按该年样本国家的人口加权;播放会保留此前的观测轨迹。它是 Gapminder 样本汇总,不能直接视为全洲居民的总体预期寿命。
library(dplyr)
library(gapminder)
library(plotly)
# Gapminder 每五年有一次记录;对样本国家按当年人口加权。
trend_groups <- c("Africa", "Americas", "Asia", "Europe")
trend_labels <- c(Africa = "非洲", Americas = "美洲",
Asia = "亚洲", Europe = "欧洲")
trend_cols <- c(Africa = "#C35A68", Americas = "#3477A4",
Asia = "#2F8D73", Europe = "#8868A8")
trend_years <- sort(unique(gapminder::gapminder$year))
trend_data <- gapminder::gapminder %>%
filter(as.character(continent) %in% trend_groups) %>%
mutate(continent = as.character(continent)) %>%
group_by(continent, year) %>%
summarise(
weighted_life = weighted.mean(lifeExp, w = pop),
n_countries = n(),
sample_pop = sum(pop),
.groups = "drop"
) %>%
mutate(
continent = factor(continent, levels = trend_groups),
hover_label = sprintf(
"%s · %d 年<br>样本人口加权预期寿命:%.1f 岁<br>样本国家或地区:%d 个<br>样本人口合计:%.1f 百万人",
trend_labels[as.character(continent)], year, weighted_life,
n_countries, sample_pop / 1e6
)
) %>%
arrange(continent, year)
# 初始状态只显示 1952 年;每一帧追加截至该年的全部观测点。
trend_plot <- plotly::plot_ly(height = 570)
for (group_name in trend_groups) {
first_point <- trend_data %>%
filter(continent == group_name, year == trend_years[1])
trend_plot <- trend_plot %>%
plotly::add_trace(
data = first_point, x = ~year, y = ~weighted_life,
text = ~hover_label, hoverinfo = "text",
type = "scatter", mode = "lines+markers",
line = list(color = unname(trend_cols[group_name]), width = 2.7),
marker = list(color = unname(trend_cols[group_name]), size = 7),
name = unname(trend_labels[group_name]),
showlegend = TRUE, inherit = FALSE
)
}
trend_plot <- trend_plot %>%
plotly::layout(
xaxis = list(title = "年份", range = c(1949, 2010),
tickvals = trend_years[c(TRUE, FALSE)],
showgrid = TRUE, zeroline = FALSE),
yaxis = list(title = "样本人口加权预期寿命(岁)", range = c(30, 84),
dtick = 10, zeroline = FALSE),
legend = list(orientation = "h", x = 0, y = 1.08,
xanchor = "left", yanchor = "bottom"),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif",
size = 11),
paper_bgcolor = "rgba(0,0,0,0)",
plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 68, r = 18, b = 110, t = 75)
)
# 显式设置帧,确保四条线都保留完整历史轨迹、固定颜色和坐标。
trend_plot <- plotly::plotly_build(trend_plot)
trend_plot$x$frames <- lapply(trend_years, function(frame_year) {
frame_data <- lapply(trend_groups, function(group_name) {
one_line <- trend_data %>%
filter(continent == group_name, year <= frame_year)
list(type = "scatter", mode = "lines+markers",
x = one_line$year, y = one_line$weighted_life,
text = one_line$hover_label)
})
list(name = as.character(frame_year), data = frame_data,
traces = as.list(seq_along(trend_groups) - 1L))
})
trend_plot$x$layout$sliders <- list(list(
active = 0, x = 0.08, y = -0.21, len = 0.89,
currentvalue = list(prefix = "年份:", font = list(size = 12)),
steps = lapply(trend_years, function(yr) {
list(label = as.character(yr), method = "animate",
args = list(list(as.character(yr)),
list(mode = "immediate", transition = list(duration = 0),
frame = list(duration = 700, redraw = TRUE))))
})
))
trend_plot$x$layout$updatemenus <- list(list(
type = "buttons", direction = "left", showactive = FALSE,
x = 1, y = 1.16, xanchor = "right", yanchor = "top",
buttons = list(
list(label = "播放", method = "animate",
args = list(NULL,
list(fromcurrent = TRUE, mode = "immediate",
transition = list(duration = 0),
frame = list(duration = 700, redraw = TRUE)))),
list(label = "暂停", method = "animate",
args = list(list(NULL),
list(mode = "immediate", transition = list(duration = 0),
frame = list(duration = 0, redraw = TRUE))))
)
))
plotly::config(trend_plot, displaylogo = FALSE, responsive = TRUE)左图在同一坐标上叠加直方图与 KDE,查看 2007 年各洲国家寿命的分布形状;右图可切换三个年份,核对中位数、四分位距与原始国家点。两图都让每个国家或地区等权,大洋洲因样本过少未纳入。
B1 · 直方图与 KDE
柱高与曲线均为概率密度;组距固定 4 岁,KDE 带宽固定 3 岁。切换大洲时,坐标范围不变,便于比较分布宽窄。
library(dplyr)
library(gapminder)
library(plotly)
hist_groups <- c("Africa", "Americas", "Asia", "Europe")
hist_labels <- c(Africa = "非洲", Americas = "美洲", Asia = "亚洲", Europe = "欧洲")
hist_cols <- c(Africa = "#C35A68", Americas = "#3477A4",
Asia = "#2F8D73", Europe = "#8868A8")
hist_data <- gapminder %>%
filter(year == 2007, as.character(continent) %in% hist_groups) %>%
mutate(continent = factor(as.character(continent), levels = hist_groups))
# 为各洲预先计算同一带宽的 KDE;切换时只改变可见图层。
hist_kde <- lapply(hist_groups, function(group_name) {
x <- hist_data$lifeExp[hist_data$continent == group_name]
density(x, bw = 3, from = 20, to = 92, n = 512)
})
names(hist_kde) <- hist_groups
hist_breaks <- seq(20, 92, by = 4)
hist_y_max <- max(vapply(hist_groups, function(group_name) {
x <- hist_data$lifeExp[hist_data$continent == group_name]
max(c(hist(x, breaks = hist_breaks, plot = FALSE)$density,
hist_kde[[group_name]]$y))
}, numeric(1))) * 1.13
# KDE 只用很浅的同色填充,保留柱形轮廓和曲线的双重线索。
rgba <- function(hex, alpha) {
rgb <- grDevices::col2rgb(hex)[, 1]
sprintf("rgba(%d,%d,%d,%.2f)", rgb[1], rgb[2], rgb[3], alpha)
}
hist_plot <- plotly::plot_ly(height = 420)
for (i in seq_along(hist_groups)) {
group_name <- hist_groups[i]
x <- hist_data$lifeExp[hist_data$continent == group_name]
dens <- hist_kde[[group_name]]
hist_plot <- hist_plot %>%
plotly::add_histogram(
x = x, histnorm = "probability density",
xbins = list(start = 20, end = 92, size = 4),
marker = list(color = hist_cols[group_name], opacity = 0.50,
line = list(color = hist_cols[group_name], width = 0.7)),
visible = i == 1, showlegend = FALSE, inherit = FALSE,
hovertemplate = "组中心:%{x:.0f} 岁<br>密度:%{y:.3f}<extra></extra>"
) %>%
plotly::add_trace(
x = dens$x, y = dens$y, type = "scatter", mode = "lines",
fill = "tozeroy", fillcolor = rgba(hist_cols[group_name], 0.10),
line = list(color = hist_cols[group_name], width = 3),
visible = i == 1, showlegend = FALSE, inherit = FALSE,
hovertemplate = "预期寿命:%{x:.1f} 岁<br>KDE:%{y:.3f}<extra></extra>"
)
}
# 每个菜单选项同时打开该洲的直方图和 KDE,其他图层保持隐藏。
hist_buttons <- lapply(seq_along(hist_groups), function(i) {
visible <- rep(FALSE, 2L * length(hist_groups))
visible[(2L * i - 1L):(2L * i)] <- TRUE
list(
label = paste0(hist_labels[hist_groups[i]], "(n = ",
sum(hist_data$continent == hist_groups[i]), ")"),
method = "update", args = list(list(visible = visible))
)
})
hist_plot %>%
plotly::layout(
xaxis = list(title = "出生时预期寿命(岁)", range = c(20, 92)),
yaxis = list(title = "国家或地区的概率密度", range = c(0, hist_y_max)),
bargap = 0.04,
updatemenus = list(list(type = "dropdown", buttons = hist_buttons,
showactive = FALSE,
x = 0.98, y = 1.15, xanchor = "right")),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif"),
paper_bgcolor = "rgba(0,0,0,0)", plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 70, r = 25, b = 55, t = 65)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)B2 · 交互箱线图与国家点
箱体显示四分位数与中位数;悬停国家点可查看原始寿命,菱形标示中位数。年份切换只改变数据,轴范围保持固定。
library(dplyr)
library(gapminder)
library(plotly)
# 同一国家在三个截面各出现一次;国家等权,不按人口加权。
box_years <- c(1952L, 1977L, 2007L)
box_groups <- c("Africa", "Americas", "Asia", "Europe")
box_labels <- c(Africa = "非洲", Americas = "美洲",
Asia = "亚洲", Europe = "欧洲")
box_colors <- c(Africa = "#C35A68", Americas = "#3477A4",
Asia = "#2F8D73", Europe = "#8868A8")
box_data <- gapminder::gapminder %>%
filter(year %in% box_years,
as.character(continent) %in% box_groups)
box_rgba <- function(hex, alpha) {
rgb <- grDevices::col2rgb(hex)[, 1]
sprintf("rgba(%d,%d,%d,%.2f)", rgb[1], rgb[2], rgb[3], alpha)
}
# 每个年份和大洲各有箱体、原始点、中位数标记三个图层。
# 原始点的纵向微偏移由排序位置计算,切换年份时位置保持一致。
box_plot <- plotly::plot_ly(height = 420)
for (year_i in seq_along(box_years)) {
this_year <- box_years[year_i]
for (group_i in seq_along(box_groups)) {
group_name <- box_groups[group_i]
group_data <- box_data %>%
filter(year == this_year, as.character(continent) == group_name) %>%
arrange(country)
values <- group_data$lifeExp
n <- length(values)
y_pos <- length(box_groups) - group_i + 1L
point_y <- y_pos + sin(seq_len(n) * 2.39996323) * 0.16
quartiles <- as.numeric(stats::quantile(
values, probs = c(0.25, 0.50, 0.75), names = FALSE
))
label <- box_labels[[group_name]]
color <- box_colors[[group_name]]
is_visible <- this_year == 2007L
stats_text <- sprintf(
"%d 年 · %s(n = %d)<br>Q1:%.1f 岁<br>中位数:%.1f 岁<br>Q3:%.1f 岁",
this_year, label, n, quartiles[1], quartiles[2], quartiles[3]
)
box_plot <- box_plot %>%
plotly::add_boxplot(
x = values, y = rep(y_pos, n), orientation = "h",
boxpoints = FALSE, quartilemethod = "linear", width = 0.46,
line = list(color = color, width = 2),
fillcolor = box_rgba(color, 0.24),
hoverinfo = "skip", visible = is_visible,
showlegend = FALSE, inherit = FALSE
) %>%
plotly::add_markers(
x = values, y = point_y,
customdata = as.character(group_data$country),
marker = list(color = box_rgba(color, 0.72), size = 5),
hovertemplate = paste0(
this_year, " 年 · ", label,
"<br>国家或地区:%{customdata}",
"<br>出生时预期寿命:%{x:.1f} 岁<extra></extra>"
),
visible = is_visible, showlegend = FALSE, inherit = FALSE
) %>%
plotly::add_markers(
x = quartiles[2], y = y_pos,
marker = list(symbol = "diamond", size = 9, color = "#FFFFFF",
line = list(color = color, width = 2)),
text = stats_text, hovertemplate = "%{text}<extra></extra>",
visible = is_visible, showlegend = FALSE, inherit = FALSE
)
}
}
# 只切换对应年份的 12 个图层,四洲的顺序与横轴范围固定。
box_buttons <- lapply(seq_along(box_years), function(i) {
visible <- rep(FALSE, 3L * length(box_groups) * length(box_years))
first_trace <- (i - 1L) * 3L * length(box_groups) + 1L
visible[first_trace:(first_trace + 3L * length(box_groups) - 1L)] <- TRUE
list(label = paste0(box_years[i], " 年"), method = "update",
args = list(list(visible = visible)))
})
box_plot <- box_plot %>%
plotly::layout(
xaxis = list(title = "出生时预期寿命(岁)", range = c(20, 92),
tickvals = c(30, 45, 60, 75, 90)),
yaxis = list(title = "", range = c(0.5, 4.5),
tickvals = 1:4,
ticktext = unname(rev(box_labels[box_groups])),
showgrid = FALSE, zeroline = FALSE),
updatemenus = list(list(
type = "dropdown", buttons = box_buttons, active = 2L,
showactive = TRUE, x = 0.98, y = 1.17, xanchor = "right"
)),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif"),
paper_bgcolor = "rgba(0,0,0,0)",
plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 67, r = 20, b = 55, t = 65)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)
box_plot左图按年份播放,观察四洲国家预期寿命在各 5
岁区间的洲内样本占比;右图用 mtcars 的 32
款车型比较五项指标的两两 Pearson 相关。右图仅保留 10
组不重复的下三角方格:格内为相关系数 \(r\),悬停或点按可查看检验结果。两图的数据与色标口径不同。
C1 · 逐年分箱占比
颜色较深的区间长期总体向较高寿命段移动,个别年份可能回退。分母是该洲该年的样本国家数;色标固定为 0—70%。
library(dplyr)
library(gapminder)
library(plotly)
# 按洲、年份与 5 岁区间统计国家或地区记录数;补齐空格。
heat_continents <- c("Africa", "Americas", "Asia", "Europe")
heat_labels <- c(Africa = "非洲", Americas = "美洲",
Asia = "亚洲", Europe = "欧洲")
age_breaks <- seq(20, 90, by = 5)
age_centres <- head(age_breaks, -1) + 2.5
heat_years <- sort(unique(gapminder::gapminder$year))
heat_source <- gapminder::gapminder %>%
filter(as.character(continent) %in% heat_continents) %>%
mutate(
continent = factor(as.character(continent), levels = heat_continents),
age_bin = cut(lifeExp, breaks = age_breaks, right = FALSE,
labels = age_centres, include.lowest = TRUE),
age_bin = as.numeric(as.character(age_bin))
)
stopifnot(!anyNA(heat_source$age_bin))
heat_counts <- heat_source %>% count(year, continent, age_bin, name = "n_bin")
heat_denominators <- heat_source %>%
count(year, continent, name = "n_continent")
heat_grid <- expand.grid(
year = heat_years, continent = heat_continents,
age_bin = age_centres, KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
) %>%
mutate(continent = factor(continent, levels = heat_continents)) %>%
left_join(heat_counts, by = c("year", "continent", "age_bin")) %>%
left_join(heat_denominators, by = c("year", "continent")) %>%
mutate(
n_bin = coalesce(n_bin, 0L),
share = n_bin / n_continent,
continent_zh = unname(heat_labels[as.character(continent)]),
lower_age = age_bin - 2.5,
upper_age = age_bin + 2.5,
hover_label = sprintf(
"%s · %s 年<br>预期寿命 %.0f—<%.0f 岁<br>%d / %d 个国家或地区(%.1f%%)",
continent_zh, year, lower_age, upper_age,
n_bin, n_continent, share * 100
)
)
# Plotly 热力图每一帧需要行列固定的二维矩阵。
heat_matrix <- function(year_value, column_name) {
one_year <- heat_grid %>%
filter(year == year_value) %>%
arrange(continent, age_bin)
matrix(one_year[[column_name]],
nrow = length(heat_continents), ncol = length(age_centres),
byrow = TRUE)
}
z_frames <- lapply(heat_years, heat_matrix, column_name = "share")
text_frames <- lapply(heat_years, heat_matrix, column_name = "hover_label")
heat_plot <- plotly::plot_ly(
x = age_centres, y = seq_along(heat_continents), z = z_frames[[1]],
type = "heatmap",
colors = c("#F1F6F7", "#C8E4E2", "#8CC8C7",
"#4A9FA9", "#40759C", "#4E567F"),
zmin = 0, zmax = 0.70, zsmooth = FALSE,
text = text_frames[[1]], hoverinfo = "text",
xgap = 1.5, ygap = 1.5,
colorbar = list(
title = "洲内样本占比",
tickvals = c(0, 0.1, 0.2, 0.3, 0.5, 0.7),
ticktext = c("0%", "10%", "20%", "30%", "50%", "70%"),
len = 0.78, thickness = 14
),
height = 500
) %>%
plotly::layout(
xaxis = list(title = "出生时预期寿命(岁;每格 5 岁)",
range = c(19.9, 90.1), tickvals = seq(20, 90, 10),
showgrid = FALSE, zeroline = FALSE),
yaxis = list(title = NULL, range = c(4.5, 0.5),
tickvals = seq_along(heat_continents),
ticktext = unname(heat_labels[heat_continents]),
showgrid = FALSE, zeroline = FALSE),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif",
size = 12),
paper_bgcolor = "rgba(0,0,0,0)", plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 67, r = 110, b = 155, t = 58)
)
# 年份帧只更新格子的占比与悬停文字,色标和坐标不变。
heat_plot <- plotly::plotly_build(heat_plot)
heat_plot$x$frames <- lapply(seq_along(heat_years), function(i) {
list(name = as.character(heat_years[i]),
data = list(list(type = "heatmap", z = z_frames[[i]],
text = text_frames[[i]])),
traces = list(0))
})
heat_plot$x$layout$sliders <- list(list(
active = 0, x = 0.05, y = -0.20, len = 0.90,
currentvalue = list(prefix = "年份:", font = list(size = 13)),
steps = lapply(heat_years, function(yr) {
list(label = as.character(yr), method = "animate",
args = list(list(as.character(yr)),
list(mode = "immediate", transition = list(duration = 0),
frame = list(duration = 800, redraw = TRUE))))
})
))
heat_plot$x$layout$updatemenus <- list(list(
type = "buttons", direction = "left", showactive = FALSE,
x = 0, y = 1.14, xanchor = "left", yanchor = "top",
buttons = list(
list(label = "播放", method = "animate",
args = list(NULL,
list(fromcurrent = TRUE, mode = "immediate",
transition = list(duration = 0),
frame = list(duration = 850, redraw = TRUE)))),
list(label = "暂停", method = "animate",
args = list(list(NULL),
list(mode = "immediate", transition = list(duration = 0),
frame = list(duration = 0, redraw = TRUE))))
)
))
plotly::config(heat_plot, displaylogo = FALSE, responsive = TRUE)C2 · 五项汽车指标相关性
Pearson r:紫色为负、绿色为正。*、**、*** 分别表示 Holm 校正后 p < 0.05、0.01、0.001;ns 为未达 0.05。悬停查看具体数值;相关不表示因果。
library(plotly)
# 1. 数据与检验:五个指标共有 32 条完整车型记录。
corr_metrics <- c("mpg", "hp", "qsec", "wt", "drat")
corr_names <- c(
mpg = "燃油经济性 mpg(英里/加仑)",
hp = "发动机功率 hp(马力)",
qsec = "四分之一英里用时 qsec(秒)",
wt = "车重 wt(千磅)",
drat = "后轴速比 drat"
)
corr_pairs <- combn(corr_metrics, 2, simplify = FALSE)
corr_tests <- lapply(corr_pairs, function(pair) {
stats::cor.test(mtcars[[pair[1]]], mtcars[[pair[2]]], method = "pearson")
})
corr_results <- data.frame(
col = vapply(corr_pairs, `[`, character(1), 1),
row = vapply(corr_pairs, `[`, character(1), 2),
r = vapply(corr_tests, function(test) unname(test$estimate), numeric(1)),
p_raw = vapply(corr_tests, function(test) test$p.value, numeric(1))
)
# 十次两两检验共用一次 Holm 多重比较校正。
corr_results$p_holm <- p.adjust(corr_results$p_raw, method = "holm")
corr_results$mark <- ifelse(corr_results$p_holm < 0.001, "***",
ifelse(corr_results$p_holm < 0.01, "**",
ifelse(corr_results$p_holm < 0.05, "*", "ns")))
# 2. 矩阵:保留下三角;上三角与对角线留空,不重复展示。
corr_z <- matrix(NA_real_, nrow = 4, ncol = 4)
corr_hover <- matrix("", nrow = 4, ncol = 4)
corr_display <- matrix("", nrow = 4, ncol = 4)
corr_p_label <- function(p_value) {
sub("^<", "<", trimws(format.pval(p_value, digits = 3, eps = 0.001)))
}
for (i in seq_len(nrow(corr_results))) {
row_id <- match(corr_results$row[i], corr_metrics) - 1L
col_id <- match(corr_results$col[i], corr_metrics)
corr_z[row_id, col_id] <- corr_results$r[i]
corr_display[row_id, col_id] <- sprintf("%.2f", corr_results$r[i])
corr_hover[row_id, col_id] <- sprintf(
"%s × %s<br>Pearson r = %+.3f<br>双侧检验 p = %s<br>Holm 校正 p = %s(%s)<br>n = 32 款车型",
corr_names[corr_results$row[i]], corr_names[corr_results$col[i]],
corr_results$r[i],
corr_p_label(corr_results$p_raw[i]),
corr_p_label(corr_results$p_holm[i]),
corr_results$mark[i]
)
}
# 3. 绘图:紫色为负相关,绿色为正相关;全图固定 -1 至 1 色标。
corr_plot <- plotly::plot_ly(
x = 1:4, y = 2:5, z = corr_z, type = "heatmap",
colorscale = list(c(0, "#B18ECA"), c(0.5, "#F2F1EE"),
c(1, "#7CB99D")),
zmin = -1, zmax = 1, zsmooth = FALSE, hoverongaps = FALSE,
text = corr_hover, hoverinfo = "text", texttemplate = corr_display,
textfont = list(color = "#233442", size = 12),
xgap = 2, ygap = 2,
colorbar = list(title = "Pearson r", thickness = 12, len = 0.66,
tickvals = c(-1, -0.5, 0, 0.5, 1)),
height = 500
) %>%
plotly::layout(
xaxis = list(title = NULL, range = c(0.5, 4.5),
tickvals = 1:4, ticktext = corr_metrics[1:4],
showgrid = FALSE, zeroline = FALSE, side = "bottom"),
yaxis = list(title = NULL, range = c(5.5, 1.5),
tickvals = 2:5, ticktext = corr_metrics[2:5],
showgrid = FALSE, zeroline = FALSE,
scaleanchor = "x", scaleratio = 1),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif",
color = "#263B4C", size = 11),
hoverlabel = list(bgcolor = "#FFFFFF", bordercolor = "#AABCC6",
font = list(color = "#263B4C", size = 11)),
paper_bgcolor = "rgba(0,0,0,0)",
plot_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 62, r = 80, b = 55, t = 20)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)
corr_plot两图都可以拖动旋转、滚轮缩放。左图用 R 自带的 volcano
高程矩阵呈现地形及等高线;右图把 2007
年四洲国家寿命的核密度曲线竖起为半透明山脊。三维位置服务于观察,各轴含义应分别读取。
D1 · 火山地形与等高线
原数据按 10 米网格记录;图中竖向比例约放大两倍,画面坡度不等于真实坡度。可悬停查看高程。
library(plotly)
# volcano 为 87 行 × 61 列的高程矩阵;坐标按 10 米网格构建。
x_grid <- seq(0, by = 10, length.out = ncol(volcano))
y_grid <- seq(0, by = 10, length.out = nrow(volcano))
plotly::plot_ly(
x = x_grid, y = y_grid, z = volcano,
type = "surface", colors = viridisLite::viridis(12),
cmin = 90, cmax = 195,
contours = list(z = list(
show = TRUE, start = 100, end = 190, size = 10,
color = "rgba(32,52,72,0.66)", width = 1.4,
usecolormap = FALSE, project = list(z = TRUE)
)),
lighting = list(ambient = 0.70, diffuse = 0.88,
roughness = 0.9, specular = 0.12, fresnel = 0.03),
lightposition = list(x = 800, y = -200, z = 1200),
colorbar = list(title = "高程(米)", len = 0.70,
thickness = 14, outlinewidth = 0),
height = 490,
hovertemplate = "东西向:%{x} 米<br>南北向:%{y} 米<br>高程:%{z} 米<extra></extra>"
) %>%
plotly::layout(
scene = list(
xaxis = list(title = "东西向距离(米)", backgroundcolor = "#F4F7F8",
showbackground = TRUE, gridcolor = "#D9E2E8"),
yaxis = list(title = "南北向距离(米)", backgroundcolor = "#F4F7F8",
showbackground = TRUE, gridcolor = "#D9E2E8"),
zaxis = list(title = "高程(米)", range = c(85, 200),
backgroundcolor = "#F4F7F8", showbackground = TRUE,
gridcolor = "#D9E2E8"),
aspectmode = "manual", aspectratio = list(x = 0.70, y = 1, z = 0.27),
camera = list(eye = list(x = 1.55, y = -1.75, z = 1.25))
),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif"),
paper_bgcolor = "rgba(0,0,0,0)",
margin = list(l = 5, r = 48, b = 5, t = 5)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)D2 · 四洲寿命三维山脊
每座山脊是一条国家等权的 KDE 曲线,统一带宽为 3.5 岁;菱形表示该洲中位数。曲线高度是核密度,山脊间的距离仅用于分组。
library(dplyr)
library(gapminder)
library(plotly)
ridge_groups <- c("Africa", "Americas", "Asia", "Europe")
ridge_names <- c(Africa = "非洲", Americas = "美洲", Asia = "亚洲", Europe = "欧洲")
ridge_cols <- c(Africa = "#C35A68", Americas = "#3477A4", Asia = "#2F8D73", Europe = "#8868A8")
ridge_source <- gapminder %>% filter(year == 2007, as.character(continent) %in% ridge_groups)
ridge_plot <- plotly::plot_ly(height = 490)
ridge_peaks <- numeric(length(ridge_groups))
for (group_id in rev(seq_along(ridge_groups))) {
group_name <- ridge_groups[group_id]
values <- ridge_source %>% filter(as.character(continent) == group_name) %>% pull(lifeExp)
curve <- stats::density(values, bw = 3.5, from = 25, to = 90, n = 80)
n <- length(curve$x)
segment <- 0:(n - 2L)
ridge_peaks[group_id] <- max(curve$y)
ridge_plot <- ridge_plot %>%
plotly::add_trace(
type = "mesh3d", x = rep(curve$x, 2L), y = rep(group_id, 2L * n),
z = c(curve$y, rep(0, n)),
i = c(segment, segment),
j = c(segment + 1L, n + segment + 1L),
k = c(n + segment + 1L, n + segment),
color = unname(ridge_cols[group_name]), opacity = 0.28,
flatshading = TRUE,
lighting = list(ambient = 0.9, diffuse = 0.2, specular = 0),
hoverinfo = "skip", showlegend = FALSE, inherit = FALSE
) %>%
plotly::add_trace(
type = "scatter3d", mode = "lines",
x = curve$x, y = rep(group_id, n), z = curve$y,
line = list(color = unname(ridge_cols[group_name]), width = 6),
name = unname(ridge_names[group_name]), showlegend = FALSE,
hovertemplate = paste0(unname(ridge_names[group_name]),
"<br>预期寿命:%{x:.1f} 岁<br>核密度:%{z:.3f}<extra></extra>"),
inherit = FALSE
)
median_value <- stats::median(values)
median_height <- stats::approx(curve$x, curve$y, xout = median_value)$y
ridge_plot <- ridge_plot %>%
plotly::add_trace(
type = "scatter3d", mode = "markers",
x = median_value, y = group_id, z = median_height,
marker = list(color = unname(ridge_cols[group_name]), size = 5,
symbol = "diamond", line = list(color = "#FFFFFF", width = 1)),
hovertemplate = paste0(unname(ridge_names[group_name]),
"<br>国家寿命中位数:%{x:.1f} 岁<extra></extra>"),
showlegend = FALSE, inherit = FALSE
)
}
ridge_plot <- ridge_plot %>%
plotly::layout(
scene = list(
xaxis = list(title = "出生时预期寿命(岁)", range = c(25, 90),
backgroundcolor = "#F4F7F8", showbackground = TRUE,
gridcolor = "#D9E2E8"),
yaxis = list(title = "大洲", range = c(0.6, 4.4),
tickvals = seq_along(ridge_groups),
ticktext = unname(ridge_names[ridge_groups]),
backgroundcolor = "#F4F7F8", showbackground = TRUE,
gridcolor = "#D9E2E8"),
zaxis = list(title = "核密度", range = c(0, max(ridge_peaks) * 1.18),
backgroundcolor = "#F4F7F8", showbackground = TRUE,
gridcolor = "#D9E2E8"),
aspectmode = "manual", aspectratio = list(x = 1.4, y = 0.82, z = 0.65),
camera = list(eye = list(x = 1.55, y = -1.8, z = 1.2))
),
font = list(family = "Arial, Microsoft YaHei, PingFang SC, sans-serif", size = 11),
paper_bgcolor = "rgba(0,0,0,0)", margin = list(l = 0, r = 0, b = 0, t = 0)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)
ridge_plot上面四组八幅图的数据与计算结果在 Knit 时已经确定。下面的 Shiny
应用则在年份、大洲或回归线开关改变后,重新筛选数据并拟合当前样本。将完整代码保存为单独文件夹内的
app.R,在 RStudio 点击 Run
App,或在该文件夹运行
shiny::runApp(".")。此应用需要运行中的 R 会话,因此课程静态
HTML 只展示代码。
library(shiny)
library(bslib)
library(plotly)
library(dplyr)
library(gapminder)
continent_cols <- c(
Africa = "#C35A68", Americas = "#3477A4", Asia = "#2F8D73",
Europe = "#8868A8", Oceania = "#B28B32"
)
continent_choices <- c(
"全部" = "all", "非洲" = "Africa", "美洲" = "Americas",
"亚洲" = "Asia", "欧洲" = "Europe", "大洋洲" = "Oceania"
)
ui <- bslib::page_sidebar(
title = "预期寿命与人均 GDP",
sidebar = bslib::sidebar(
shiny::sliderInput("selected_year", "年份", min = 1952, max = 2007,
value = 2007, step = 5,
animate = shiny::animationOptions(interval = 1100)),
shiny::selectInput("selected_continent", "大洲",
choices = continent_choices, selected = "all"),
shiny::checkboxInput("show_fit", "显示当前样本回归线", value = TRUE)
),
bslib::card(
plotly::plotlyOutput("scatter_plot", height = "510px"),
shiny::textOutput("sample_note")
),
theme = bslib::bs_theme(version = 5, bootswatch = "flatly")
)
server <- function(input, output, session) {
selected_data <- shiny::reactive({
data <- gapminder %>% filter(year == input$selected_year)
if (input$selected_continent != "all") {
data <- data %>% filter(as.character(continent) == input$selected_continent)
}
data
})
output$sample_note <- shiny::renderText({
data <- selected_data()
paste0(
"当前显示 ", nrow(data), " 个国家或地区。",
if (nrow(data) < 3L) "样本少于 3 个,不绘制回归线。" else
if (isTRUE(input$show_fit))
"回归线仅概括当前样本的关联,不能解释为因果效应。" else
"回归线已隐藏,可用左侧开关重新显示。"
)
})
output$scatter_plot <- plotly::renderPlotly({
data <- selected_data()
data$continent <- factor(as.character(data$continent), levels = names(continent_cols))
plot <- plotly::plot_ly(
data, x = ~gdpPercap, y = ~lifeExp,
color = ~continent, colors = unname(continent_cols),
type = "scatter", mode = "markers",
marker = list(size = 9, opacity = 0.68),
text = ~paste0(
country, "<br>", year, " 年<br>人均 GDP:",
format(round(gdpPercap), big.mark = ",", scientific = FALSE),
" 美元<br>预期寿命:", sprintf("%.1f", lifeExp),
" 岁<br>人口:", sprintf("%.1f", pop / 1e6), " 百万人"
), hoverinfo = "text"
)
if (isTRUE(input$show_fit) && nrow(data) >= 3L) {
fit <- lm(lifeExp ~ log10(gdpPercap), data = data)
x_grid <- exp(seq(log(min(data$gdpPercap)),
log(max(data$gdpPercap)), length.out = 80))
fit_data <- data.frame(
gdpPercap = x_grid,
lifeExp = as.numeric(predict(fit, data.frame(gdpPercap = x_grid)))
)
plot <- plot %>% plotly::add_lines(
data = fit_data, x = ~gdpPercap, y = ~lifeExp,
inherit = FALSE, line = list(color = "#263B4C", width = 2.5),
name = "当前样本回归线", hoverinfo = "skip"
)
}
plot %>%
plotly::layout(
xaxis = list(title = "人均 GDP(美元,对数刻度)", type = "log",
range = c(2.35, 5.1)),
yaxis = list(title = "出生时预期寿命(岁)", range = c(20, 88)),
legend = list(orientation = "v"),
margin = list(l = 65, r = 30, b = 55, t = 20)
) %>%
plotly::config(displaylogo = FALSE, responsive = TRUE)
})
}
shiny::shinyApp(ui, server)bslib
用于应用的响应式页面布局。当只选大洋洲时,每年仅有两个国家,应用仍显示散点,但不会拟合回归线。
读懂应用的三个位置:
sliderInput()、selectInput() 和
checkboxInput() 接收操作,数值保存在 input
中。reactive()
根据这些输入重新筛选记录;选项改变后,依赖它的计算会更新。renderPlotly() 重新生成图,plotlyOutput()
指定图在页面上的位置。进一步阅读: Plotly for R · Shiny 官方入门 · shinylive for R · gganimate。
同一批对象被多次测量时,动画可以把各次观测连起来,帮助跟踪变化。本例用柱长呈现体重增长,用位置交换呈现排名变化;横轴固定,各帧可按同一尺度比较。
R 自带的 ChickWeight 记录 50 只鸡雏在四种饲料下的体重,每行是一只鸡雏在某个日龄的测量。本例保留 0—21 天内 12 次记录齐全的 45 只,展示体重前十名,并用全部 45 只计算均值。
weight 体重(克)Time 日龄(天)Chick 鸡雏编号Diet 饲料类型测量日之间采用线性过渡,过渡帧不代表新增实测记录。同重时按编号排列。本图用于观察个体变化,不能仅凭前十名判断哪种饲料效果更好。
首次安装:在 Console 中运行一次。
包调用:写入作图脚本,每次运行时执行。
gifski 合成 GIF,ragg
渲染画面,systemfonts 检查字体;完整脚本通过
包名::函数名()
调用这三个辅助包。生成动画需要等待渲染完成,时间取决于帧数与图片尺寸。
代码中的 group = chick_key
保持同一只鸡雏在各帧中的对应关系;transition_time(age_day)
按实际日龄推进,ease_aes("linear") 控制过渡方式。最后由
animate() 渲染、anim_save() 保存 GIF。JAMA
色盘的四个 HEX 写在 palette_colors
中,palette_alpha = 0.84
统一控制柱体和图例的透明度。字体可通过 font_family
改为本机已有字体。
# EasyPlot 动画测试:R 内置 ChickWeight 的鸡雏体重排名。
# 参考公众号的横向柱状竞速思路,重新实现数据处理、固定坐标和饲料配色。
# 一次性依赖配置可在 R Console 运行:
# install.packages(c("ggplot2", "dplyr", "gganimate", "gifski", "ragg", "systemfonts"))
# 在终端运行:Rscript build-chickweight-race.R
if (.Platform$OS.type == "windows") Sys.setlocale("LC_CTYPE", ".UTF-8")
task_library <- Sys.getenv("EASYPLOT_ANIMATION_LIB", unset = "")
if (nzchar(task_library)) .libPaths(c(task_library, .libPaths()))
required_packages <- c("ggplot2", "dplyr", "gganimate", "gifski", "ragg", "systemfonts")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages)) {
stop("请先安装缺少的包:", paste(missing_packages, collapse = ", "))
}
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(gganimate)
})
script_argument <- grep("^--file=", commandArgs(), value = TRUE)
script_dir <- if (length(script_argument)) {
dirname(normalizePath(sub("^--file=", "", script_argument[1]), winslash = "/"))
} else getwd()
output_dir <- file.path(script_dir, "output")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# 1. 保留真实测量。Chick 是有序因子,先转字符再转整数才是鸡雏编号。
raw_data <- datasets::ChickWeight %>%
transmute(
chick_id = as.integer(as.character(Chick)),
age_day = Time,
weight_g = weight,
diet = as.character(Diet)
)
measurement_days <- sort(unique(raw_data$age_day))
stopifnot(nrow(raw_data) == 578L, !anyNA(raw_data),
!anyDuplicated(raw_data[c("chick_id", "age_day")]))
# 2. 本样板比较同一批鸡雏。45 只具有全部 12 次记录;5 只后续记录不全。
# 退出原因未提供,不能将记录不全解释为死亡;这些数据不作外推或补零。
complete_ids <- raw_data %>%
count(chick_id, name = "observations") %>%
filter(observations == length(measurement_days)) %>%
pull(chick_id)
excluded_ids <- setdiff(unique(raw_data$chick_id), complete_ids)
complete_data <- raw_data %>% filter(chick_id %in% complete_ids)
rank_data <- complete_data %>%
group_by(age_day) %>%
arrange(desc(weight_g), chick_id, .by_group = TRUE) %>%
mutate(rank = row_number()) %>%
ungroup() %>%
mutate(
chick_key = factor(chick_id),
chick_label = paste0("鸡雏 ", chick_id),
diet_label = factor(diet, levels = as.character(1:4),
labels = paste0("饲料 ", 1:4)),
display_y = -rank
)
stopifnot(length(complete_ids) == 45L, nrow(rank_data) == 540L,
all(table(rank_data$age_day) == 45L))
# 同重时按鸡雏编号升序排位置,并不表示体重存在差异。
# 先在全部 45 只的相邻实测之间插值,每 0.2 天重新排名,再取前 10。
# 这样避免只对测量日已入榜的鸡雏做过渡,漏掉测量间进入前十的个体。
animation_days <- sort(unique(c(round(seq(0, 21, by = 0.2), 6), measurement_days)))
interpolated_data <- bind_rows(lapply(split(complete_data, complete_data$chick_id),
function(chick_data) {
data.frame(
chick_id = chick_data$chick_id[1],
age_day = animation_days,
weight_g = approx(chick_data$age_day, chick_data$weight_g,
xout = animation_days, rule = 1)$y,
diet = chick_data$diet[1]
)
}
))
animation_data <- interpolated_data %>%
group_by(age_day) %>%
arrange(desc(weight_g), chick_id, .by_group = TRUE) %>%
mutate(rank = row_number()) %>%
slice_head(n = 10L) %>%
ungroup() %>%
mutate(
chick_key = factor(chick_id, levels = sort(complete_ids)),
chick_label = paste0("鸡雏 ", chick_id),
diet_label = factor(diet, levels = as.character(1:4),
labels = paste0("饲料 ", 1:4)),
display_y = -rank
)
stopifnot(!anyNA(animation_data), all(table(animation_data$age_day) == 10L))
endpoint_check <- animation_data %>%
filter(age_day %in% measurement_days) %>%
select(age_day, chick_id, rank, weight_g) %>%
arrange(age_day, rank)
observed_check <- rank_data %>%
filter(rank <= 10L) %>%
select(age_day, chick_id, rank, weight_g) %>%
arrange(age_day, rank)
stopifnot(identical(endpoint_check, observed_check))
# 均值来自全部 45 只完整记录鸡雏,不能对展示的前 10 名取平均。
observed_means <- complete_data %>%
group_by(age_day) %>%
summarise(mean_weight_g = mean(weight_g), n_chicks = n(), .groups = "drop")
mean_data <- interpolated_data %>%
group_by(age_day) %>%
summarise(mean_weight_g = mean(weight_g), n_chicks = n(), .groups = "drop") %>%
mutate(mean_key = factor("cohort_mean"))
mean_check <- mean_data %>%
filter(age_day %in% measurement_days) %>%
select(age_day, mean_weight_g, n_chicks)
stopifnot(all(mean_data$n_chicks == 45L),
isTRUE(all.equal(mean_check, observed_means, tolerance = 1e-12)))
# EasyPlot 色盘库 ggsci.jama.default 的前四个原始色;颜色固定对应四种饲料。
# alpha 只调整柱体和图例的透明程度,不表示样本量或其他数据含义。
palette_id <- "ggsci.jama.default"
palette_colors <- c("饲料 1" = "#374E55", "饲料 2" = "#DF8F44",
"饲料 3" = "#00A1D5", "饲料 4" = "#B24745")
palette_alpha <- 0.84
font_family <- Sys.getenv("EASYPLOT_FONT", unset = "Microsoft YaHei")
font_match <- systemfonts::match_fonts(font_family)
stopifnot(file.exists(font_match$path))
if (.Platform$OS.type == "windows") {
# gganimate 的排版测量可能调用 Windows 设备,需在当前进程注册同一字体。
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(font_family)), font_family))
}
ink <- "#283A46"
# 3. 在实际绘帧时才将数值四舍五入,避免数字标签在两次测量间跳变。
# 两根柱子交换位置时,临时淡出挤在一起的标签,避免文字互相盖住。
fade_crossing_labels <- function(data) {
if (nrow(data) < 2L) return(data)
original_alpha <- ifelse(is.na(data$alpha), 1, data$alpha)
distances <- abs(outer(data$y, data$y, "-"))
diag(distances) <- Inf
distances[, original_alpha < 0.10] <- Inf
clearance <- apply(distances, 1, min)
label_alpha <- pmin(1, pmax(0, (clearance - 0.38) / 0.34))
data$alpha <- original_alpha * label_alpha
data
}
geom_race_text <- function(mapping = NULL, data = NULL, format_weight = FALSE,
format_mean = FALSE, inherit.aes = TRUE,
show.legend = NA, ...) {
race_geom <- ggproto("GeomRaceText", GeomText,
draw_panel = function(data, panel_params, coord, ...) {
data <- fade_crossing_labels(data)
if (format_weight) data$label <- paste0(round(data$label), " g")
if (format_mean) data$label <- sprintf("平均体重 %.1f g", data$label)
GeomText$draw_panel(data, panel_params, coord, ...)
}
)
layer(mapping = mapping, data = data, geom = race_geom,
stat = "identity", position = "identity", inherit.aes = inherit.aes,
show.legend = show.legend,
params = list(...))
}
make_plot <- function(plot_data, age_label, mean_plot_data) {
ggplot(plot_data,
aes(x = weight_g, y = display_y, group = chick_key)) +
geom_vline(xintercept = seq(0, 400, 100), colour = "#E4E9ED",
linewidth = 0.35) +
geom_col(aes(fill = diet_label), orientation = "y", width = 0.70,
alpha = palette_alpha, colour = "#53616A", linewidth = 0.18) +
# 白色细底线让均值标记在深浅柱子上都能看清;深色虚线才是参考线。
geom_vline(data = mean_plot_data,
aes(xintercept = mean_weight_g, group = mean_key),
colour = "white", linewidth = 1.0, show.legend = FALSE) +
geom_vline(data = mean_plot_data,
aes(xintercept = mean_weight_g, group = mean_key),
colour = ink, linetype = "dashed", linewidth = 0.48,
show.legend = FALSE) +
geom_race_text(data = mean_plot_data,
aes(x = mean_weight_g + 4, y = 0.15,
label = mean_weight_g, group = mean_key),
format_mean = TRUE, inherit.aes = FALSE,
hjust = 0, family = font_family, fontface = "bold",
size = 4.0, colour = ink, show.legend = FALSE) +
geom_race_text(aes(x = -5, label = chick_label), hjust = 1,
family = font_family, fontface = "bold", size = 4.2,
colour = ink) +
geom_race_text(aes(label = weight_g), format_weight = TRUE, hjust = -0.15,
family = font_family, fontface = "bold", size = 4.1,
colour = ink) +
scale_fill_manual(values = palette_colors, drop = FALSE,
name = "饲料类型",
guide = guide_legend(nrow = 1,
override.aes = list(alpha = palette_alpha))) +
scale_x_continuous(breaks = seq(0, 400, 100), expand = expansion(mult = 0)) +
scale_y_continuous(breaks = NULL, expand = expansion(mult = 0)) +
coord_cartesian(xlim = c(-78, 418), ylim = c(-10.58, 0.60), clip = "on") +
labs(title = age_label,
subtitle = "45 只全程测量鸡雏中,体重前 10 名",
x = "体重(g)", y = NULL,
caption = paste0(
"数据来源:R datasets::ChickWeight;Hand & Crowder (1996), Table A.2。\n",
"虚线为 45 只鸡雏的平均体重;测量间为线性过渡。")) +
theme_minimal(base_size = 14, base_family = font_family) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(colour = "#596B77", size = 11),
axis.title.x = element_text(colour = ink, size = 12, margin = margin(t = 10)),
plot.title = element_text(colour = ink, face = "bold", size = 22),
plot.subtitle = element_text(colour = "#657580", size = 12, margin = margin(b = 7)),
plot.caption = element_text(colour = "#596B77", size = 9.2,
hjust = 0, lineheight = 1.35, margin = margin(t = 12)),
plot.caption.position = "plot",
legend.position = "top", legend.justification = "left",
legend.title = element_text(colour = ink, size = 11),
legend.text = element_text(colour = ink, size = 11),
legend.key.size = grid::unit(0.40, "cm"),
legend.margin = margin(0, 0, 5, 0),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(18, 23, 14, 15)
)
}
# 4. 动画按真实日龄推进,20→21 天用时是其他两天间隔的一半。
# x(体重)、y(展示位置)和数值标签均在相邻测量之间线性过渡。
# 非测量时刻是视觉插值,不应作为新的实测数据使用。
animated_plot <- make_plot(animation_data, "日龄 {sprintf('%.1f', frame_time)} 天", mean_data) +
transition_time(age_day) +
ease_aes("linear") +
enter_fade() + exit_fade()
animation <- animate(
animated_plot, nframes = 252, fps = 18,
width = 1100, height = 800, res = 120,
device = "ragg_png", renderer = gifski_renderer(loop = TRUE),
start_pause = 18, end_pause = 36
)
animation_path <- file.path(output_dir, "chickweight_top10_race_mean.gif")
anim_save(animation_path, animation = animation)
# 5. 留下真实观测时刻的静态图、排名表和来源说明,便于核对动画。
for (age in c(0, 10, 21)) {
still_plot <- make_plot(filter(rank_data, age_day == age, rank <= 10L), paste0("日龄 ", age, " 天"),
filter(mean_data, age_day == age))
ggsave(file.path(output_dir, sprintf("chickweight_day_%02d.png", age)),
still_plot, device = ragg::agg_png, width = 1100, height = 800,
units = "px", dpi = 120, bg = "white")
}
write.csv(raw_data, file.path(output_dir, "chickweight_source_data.csv"),
row.names = FALSE, fileEncoding = "UTF-8")
write.csv(rank_data, file.path(output_dir, "chickweight_observed_rankings.csv"),
row.names = FALSE, fileEncoding = "UTF-8")
write.csv(animation_data, file.path(output_dir, "chickweight_interpolated_top10.csv"),
row.names = FALSE, fileEncoding = "UTF-8")
write.csv(observed_means, file.path(output_dir, "chickweight_observed_means.csv"),
row.names = FALSE, fileEncoding = "UTF-8")
notes <- c(
"图形:鸡雏体重前十名横向柱状竞速动画。",
"数据:R datasets::ChickWeight;体重单位为克,Time 为日龄。",
paste0("保留:", length(complete_ids), " 只具有全部 12 次测量的鸡雏,共 540 条记录。"),
paste0("记录不全而未纳入本图的编号:", paste(sort(excluded_ids), collapse = ", "), "。原因未提供。"),
"同重时按鸡雏编号升序排列。颜色固定表示饲料类型,不随排名改变。",
paste0("色盘:", palette_id, ",取前四个原始色:", paste(palette_colors, collapse = ", "), "。"),
paste0("柱体与图例透明度 alpha = ", palette_alpha,
",在白色底上保留四种饲料的固定色彩对应;透明度不编码数据。"),
"移动虚线和动态读数表示全部 45 只完整记录鸡雏的平均体重,不是前十名的均值。",
"第 0、10、21 天平均体重分别为 41.1、110.1、218.7 克。",
"横轴固定为 0—400 克;图框左部仅为标签留白。",
"非测量时刻的体重与展示位置为线性插值,仅用于动画过渡;不作新的实测数据或推断。",
"此图用于看个体体重和排名的变化,不能据前十名判断哪种饲料效果更好。",
"R 官方数据说明:https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/ChickWeight.html",
"动画思路参考:https://mp.weixin.qq.com/s/1_noj04fAWQg8_3EvFzv3g",
"原文作者:Leo Duan;公众号:AI与R语言实践。原文数据是年度通胀率,本脚本未执行原文代码。",
paste0("字体:", font_family, ";解析路径:", font_match$path),
paste0("软件:", paste(required_packages,
vapply(required_packages, function(x) as.character(packageVersion(x)), character(1)),
collapse = "; "))
)
writeLines(enc2utf8(notes), file.path(output_dir, "README.txt"), useBytes = TRUE)
writeLines(capture.output(sessionInfo()), file.path(output_dir, "session_info.txt"))
stopifnot(file.exists(animation_path), file.info(animation_path)$size > 100000)
cat("Saved animation:", normalizePath(animation_path, winslash = "/"), "\n")
cat("Complete chicks:", length(complete_ids), "; observed rows:", nrow(rank_data), "\n")数据说明:R 官方 ChickWeight 文档 · 动画思路参考:Leo Duan · AI与R语言实践
EasyPlot 是供 AI 智能体使用的科研作图 Skill,提供图形模板、配色和导出参考。默认使用 R,也可指定 Python。作者及公众号:川柏同学
先检查变量与缺失值,再选图,生成可修改的代码和图片。
先检查标注、字体、配色和导出尺寸;要修改图形,还需提供数据或源代码。
把这句话发给支持 Agent Skills 的智能体:
帮我安装这个 Skill:https://github.com/Rimagination/easyplot
在支持斜杠菜单的工具中输入 /,选择 EasyPlot,再说出数据和想画的图。
请用 EasyPlot,把 R 自带的 Titanic 数据画成流向图。
成图后核对数据汇总和图中人数;这张图的三列是并列属性,不表示时间先后。
资料:EasyPlot 项目仓库 · 公众号原文
1 · 图形与适用问题
散点图(Scatter
Plot)用一个点表示一条观测,横、纵坐标分别对应两个数值变量。手头有成对的数值数据,想看它们是否同向变化、是否分成若干群,或是否有远离其他观测的点时,可以先画散点图。例如,每个国家都有“人均
GDP”和“预期寿命”两个数值,就可以把国家放到二维坐标中。绘点函数是
geom_point()。
2 · 数据处理
先保留 gapminder 的 2007
年记录。这样每个国家或地区在图中只有一个点,避免把不同年份的观测混在同一幅横截面图里。以下数据对象也供后面的最小出图代码使用。
3 · 最小代码与基础图
把人均 GDP 放在横轴、预期寿命放在纵轴,再用 geom_point()
画点。先看基本点云,其他设置留到进阶图。
2007 年各国人均 GDP 与预期寿命的基础散点图。每个点代表一个样本国家。
读图方法
gdpPercap
是人均 GDP,纵轴 lifeExp 是出生时预期寿命(岁);每个点代表一个国家或地区。从图中可直接读出的主要信息
4 · 进阶成图与美化思路
进阶图继续使用 gapminder 的国家记录,把观察年份扩展到
1952、1977 和 2007
年,并按非洲、美洲、亚洲、欧洲分成四个面板。每个面板分别显示三个年份的散点及线性回归线:同一年份的点、线和方程使用同一种颜色,线型提供第二种辨认线索。气泡大小表示人口,半透明散点帮助辨认重叠区域。
各条回归线只覆盖相应大洲、年份实际观测到的人均 GDP
范围;同色浅带表示拟合均值的逐点 95% 置信区间。方程中的
x 是 log10(人均 GDP),y
是出生时预期寿命。拟合时每个国家或地区计一次,气泡所示人口不作为回归权重。大洋洲每年只有两个国家,未纳入本图;图中也不选择性标注个别国家。
1952、1977 与 2007 年四大洲的人均 GDP、出生时预期寿命与人口气泡图。各面板按年份绘制散点、线性回归线及拟合均值的逐点 95% 置信区间。大洋洲因每年仅有两个国家,未纳入本图。
美化思路
log10 转换后排列。plot_scatter_by_continent <- local({
library(tidyverse)
library(gapminder)
# 1. 选择本机已有的中文无衬线字体,不下载字体
os_name <- Sys.info()[["sysname"]]
font_candidates <- if (identical(os_name, "Windows")) {
c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei")
} else if (identical(os_name, "Darwin")) {
c("PingFang SC", "Heiti SC", "Hiragino Sans GB")
} else {
c("Noto Sans CJK SC", "WenQuanYi Micro Hei")
}
available_fonts <- if (requireNamespace("systemfonts", quietly = TRUE)) {
unique(systemfonts::system_fonts()[["family"]])
} else {
character()
}
matched_fonts <- font_candidates[font_candidates %in% available_fonts]
plot_font <- if (length(matched_fonts) > 0) matched_fonts[[1]] else "sans"
# 兼容 Windows 的基础绘图设备
if (.Platform$OS.type == "windows" && plot_font != "sans") {
do.call(
grDevices::windowsFonts,
setNames(
list(grDevices::windowsFont(plot_font)),
plot_font
)
)
}
# 2. 数据:三个年份、四大洲
scatter_main <- gapminder %>%
filter(
year %in% c(1952, 1977, 2007),
continent %in% c("Africa", "Americas", "Asia", "Europe")
) %>%
mutate(
year = factor(year, levels = c(1952, 1977, 2007)),
continent = factor(
as.character(continent),
levels = c("Africa", "Americas", "Asia", "Europe")
),
log_gdp = log10(gdpPercap),
pop_million = pop / 1e6
)
# 3. 每个大洲、每个年份分别拟合线性回归
fit_grid <- scatter_main %>%
group_by(continent, year) %>%
group_modify(~ {
model <- lm(lifeExp ~ log_gdp, data = .x)
new_data <- tibble(
log_gdp = seq(
min(.x$log_gdp),
max(.x$log_gdp),
length.out = 100
)
)
prediction <- predict(
model,
newdata = new_data,
interval = "confidence",
level = 0.95
)
new_data %>%
mutate(
fitted = prediction[, "fit"],
lower = prediction[, "lwr"],
upper = prediction[, "upr"],
intercept = unname(coef(model)[1]),
slope = unname(coef(model)[2])
)
}) %>%
ungroup()
# 4. 方程放在各分面的空白区域
regression_equations <- fit_grid %>%
distinct(continent, year, intercept, slope) %>%
mutate(
x_eq = if_else(continent == "Africa", 4.98, 2.48),
hjust_eq = if_else(continent == "Africa", 1, 0),
y_eq = case_when(
year == "1952" ~ 87.5,
year == "1977" ~ 84.5,
year == "2007" ~ 81.5
),
equation = sprintf(
"%s: y = %.1fx %s %.1f",
as.character(year),
slope,
if_else(intercept < 0, "-", "+"),
abs(intercept)
)
)
# 略深的同色系让细线和方程文字在白底上更清楚
year_cols <- c(
"1952" = "#A45262",
"1977" = "#3B72A1",
"2007" = "#327858"
)
year_linetypes <- c(
"1952" = "solid",
"1977" = "22",
"2007" = "42"
)
continent_labs <- c(
Africa = "非洲",
Americas = "美洲",
Asia = "亚洲",
Europe = "欧洲"
)
# 5. 成图
ggplot(scatter_main, aes(x = log_gdp, y = lifeExp)) +
# 置信带:不单独生成图例
geom_ribbon(
data = fit_grid,
aes(
x = log_gdp,
ymin = lower,
ymax = upper,
fill = year,
group = year
),
inherit.aes = FALSE,
alpha = 0.09,
colour = NA,
show.legend = FALSE
) +
# 回归线:同一年份同时对应固定颜色与线型
geom_line(
data = fit_grid,
aes(
x = log_gdp,
y = fitted,
colour = year,
linetype = year,
group = year
),
inherit.aes = FALSE,
linewidth = 0.8
) +
# 散点:人口控制面积,透明度帮助辨认重叠
geom_point(
aes(size = pop_million, colour = year),
shape = 16,
alpha = 0.58
) +
# 方程:与相应年份的回归线同色
geom_text(
data = regression_equations,
aes(
x = x_eq,
y = y_eq,
label = equation,
colour = year,
hjust = hjust_eq
),
inherit.aes = FALSE,
family = plot_font,
fontface = "bold",
size = 3,
show.legend = FALSE
) +
facet_wrap(
~ continent,
ncol = 2,
labeller = as_labeller(continent_labs)
) +
scale_colour_manual(
name = "年份",
values = year_cols,
breaks = names(year_cols)
) +
scale_linetype_manual(
name = "年份",
values = year_linetypes,
breaks = names(year_linetypes)
) +
scale_fill_manual(
values = year_cols,
guide = "none"
) +
scale_size(
name = "人口(百万人)",
range = c(1.9, 8.0),
breaks = c(100, 500, 1000)
) +
scale_x_continuous(
breaks = c(3, 4, 5),
labels = c("1,000", "10,000", "100,000")
) +
scale_y_continuous(
breaks = c(30, 45, 60, 75, 90)
) +
coord_cartesian(
xlim = c(2.4, 5.1),
ylim = c(23, 90),
clip = "on"
) +
labs(
x = "人均 GDP(对数刻度)",
y = "出生时预期寿命(岁)",
caption = paste0(
"注:方程中 x = log10(人均 GDP),y = 出生时预期寿命;",
"系数显示为一位小数。\n",
"同色浅带为拟合均值的逐点 95% 置信区间;",
"大洋洲每年仅有两个国家,未纳入分面回归。"
)
) +
# 只放大“年份”图例中的点;颜色与线型仍合并为一个图例
guides(
colour = guide_legend(
order = 1,
override.aes = list(size = 3.6, alpha = 1)
),
linetype = guide_legend(order = 1),
size = guide_legend(
order = 2,
override.aes = list(
colour = "#6B8791",
alpha = 0.75
)
)
) +
theme_classic(base_family = plot_font, base_size = 11) +
theme(
panel.grid.major.y = element_line(
colour = "#E8EEF1",
linewidth = 0.22
),
axis.line = element_line(
colour = "#738693",
linewidth = 0.3
),
axis.text = element_text(
colour = "#405362",
size = 9
),
axis.title = element_text(
colour = "#263B4C",
face = "bold",
size = 10.5
),
strip.background = element_blank(),
strip.text = element_text(
colour = "#263B4C",
face = "bold",
size = 11,
margin = margin(b = 5)
),
panel.spacing = grid::unit(0.7, "lines"),
legend.position = "right",
legend.box = "vertical",
legend.title.position = "top",
legend.title = element_text(
colour = "#263B4C",
face = "bold",
size = 9.5
),
legend.text = element_text(
colour = "#405362",
size = 9
),
legend.key.width = grid::unit(1.4, "cm"),
legend.key.height = grid::unit(0.85, "cm"),
plot.caption = element_text(
colour = "#5E6F7A",
size = 8,
lineheight = 1.35,
hjust = 0
),
plot.margin = margin(8, 12, 10, 8)
)
})
print(plot_scatter_by_continent)1 · 图形与适用问题
折线图(Line Chart)用于呈现指标随时间或有序序列变化的趋势。它通过连接相邻观测点,强调变化的方向与节奏,适合表达长期上升/下降、阶段性波动与结构性拐点等动态特征。
适合回答的问题
趋势刻画:
识别指标随时间演化的总体方向与增长/下降速率。
(如:过去数十年预期寿命的变化轨迹)
多组对比:
在同一时间轴上比较不同对象的变化幅度、速度与阶段差异。
(如:不同国家的增长路径是否同步)
拐点识别:
定位趋势发生明显转折或出现异常波动的时间点/阶段。
(Turning points / shocks)
绘图函数与分组
在 ggplot2 中,折线图对应的几何对象函数为:
geom_line()
关键逻辑:当图中包含多条线时,需要明确“哪些点属于同一条线”。最常见做法是将分组变量映射到
colour(同时完成分组);若不希望上色,也可显式设置
aes(group = ...)。
2 · 数据处理
本节使用 gapminder 数据,选取 China、United
States、Japan、India 四个国家,绘制 1952—2007
年预期寿命的变化轨迹,用于对比不同国家的长期趋势与阶段性差异。
数据与最小图的分工:先筛选四国全部观测年份,再用最少的绘图语句连接同一国家的相邻观测。color = country
同时指定分组,避免把不同国家的点接在一起。
library(tidyverse)
library(gapminder)
data_line <- gapminder %>%
filter(country %in% c("China", "United States", "Japan", "India"))3 · 最小代码与基础图
横轴为年份,纵轴为出生时预期寿命;geom_line()
画出最基本的四国趋势。
这段代码的关键映射:
year 决定点在时间轴上的位置,lifeExp
决定线的高度。color = country
同时把四个国家分成四条线;这一步若省略,连线可能把不同国家的数据接在一起。
1952—2007 年四个样本国家的预期寿命变化。每条线连接同一国家相邻年份的观测。
怎么读
据此读出的主要发现
重要
折线图把相邻观测连接起来,读者通常会据此理解变化方向。使用前应确认横轴有明确顺序,并说明连线表示的是同一对象随时间变化,还是不同群体沿有序区间的差异。
1. 无序类别变量(Nominal Data)
2. 离散区间 / 队列(Discrete Cohorts)
[10–19], [20–29],
[30–39])时,可以用连线展示各年龄段的平均值走势。使用原则:
只有当相邻 x
值的顺序与连线含义都清楚时,才适合用折线图;否则优先采用柱状图、点图或箱线图等离散比较图形。
4 · 进阶成图与美化思路
基础图用四个国家说明折线如何连接按年份排列的观测值。进阶图改以大洲为分析单位:主图比较非洲、美洲、亚洲、欧洲的长期趋势,并加入142 国样本参照线;浅色带显示同一大洲各国预期寿命的中间 50% 分布。下方四个小图则显示各洲与同期样本参照的差值。
这张图同时呈现水平、洲内差异和相对差距。完整代码收在下方,按数据汇总、配色与字体、主图和差值图分段注释。组合上下两部分使用
gridExtra;如尚未安装,可先在 Console 运行
install.packages("gridExtra")。
1952—2007 年四洲国家预期寿命的人口加权均值、洲内国家分布,以及相对 142 国样本参照的同期差值(岁)。
读图与设计要点:主图的实线是按各国总人口加权的洲内国家预期寿命均值;四洲分别用颜色和点形标识。深灰虚线与空心圆标出 Gapminder 142 国样本的人口加权参照值,包含未单独绘制的大洋洲。浅色带是当年洲内各国预期寿命的第 25—75 百分位范围,每个国家等权,表示国家间差异,不是置信区间。加权均值可能落在等权分位带之外,例如 1952 年的美洲,这是两种权重口径不同所致。下方每个小图计算“该洲加权均值 − 同年样本参照值”,单位为岁;虚线 0 表示两者相同,正负值分别表示高于或低于参照。数据每五年记录一次,线段连接观测年份,不代表中间年份已被逐年测量;样本参照也不等同于官方公布的全球预期寿命。
# 进阶折线图:四洲趋势、洲内差异与 142 国样本参照
library(dplyr)
library(ggplot2)
library(gapminder)
library(gridExtra)
# 1. 配色:Okabe-Ito 色相的深色改编,辅以点形
continent_order <- c("Africa", "Americas", "Asia", "Europe")
series_order <- c(continent_order, "World")
series_zh <- c(
Africa = "非洲", Americas = "美洲",
Asia = "亚洲", Europe = "欧洲", World = "142 国样本"
)
continent_cols <- c(
Africa = "#B45700", Americas = "#0072B2",
Asia = "#602778", Europe = "#006B4E"
)
world_col <- "#2E353B"
series_cols <- c(continent_cols, World = world_col)
series_linetypes <- c(
Africa = "solid", Americas = "solid",
Asia = "solid", Europe = "solid", World = "22"
)
series_shapes <- c(
Africa = 16, Americas = 17,
Asia = 15, Europe = 18, World = 1
)
# 2. 字体:与进阶气泡图同一选择顺序
os_name <- Sys.info()[["sysname"]]
font_candidates <- if (identical(os_name, "Windows")) {
c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei")
} else if (identical(os_name, "Darwin")) {
c("PingFang SC", "Heiti SC", "Hiragino Sans GB")
} else {
c("Noto Sans CJK SC", "WenQuanYi Micro Hei")
}
available_fonts <- if (requireNamespace("systemfonts", quietly = TRUE)) {
unique(systemfonts::system_fonts()[["family"]])
} else character()
matched_fonts <- font_candidates[font_candidates %in% available_fonts]
plot_font <- if (length(matched_fonts)) matched_fonts[[1]] else "sans"
if (.Platform$OS.type == "windows" && plot_font != "sans") {
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(plot_font)), plot_font))
}
# 3. 数据:加权均值、等权国家四分位、142 国样本基准
course_data <- gapminder::gapminder
cont_year <- course_data %>%
filter(continent %in% continent_order) %>%
group_by(continent, year) %>%
summarise(
n_country = n(),
life_w = weighted.mean(lifeExp, pop),
q25 = unname(quantile(lifeExp, 0.25)),
q75 = unname(quantile(lifeExp, 0.75)),
.groups = "drop"
) %>%
mutate(continent = factor(as.character(continent),
levels = continent_order))
world_year <- course_data %>%
group_by(year) %>%
summarise(
n_country = n(),
world_w = weighted.mean(lifeExp, pop),
.groups = "drop"
)
stopifnot(
nrow(cont_year) == 48L,
nrow(world_year) == 12L,
all(world_year$n_country == 142L),
all(cont_year$q25 <= cont_year$q75),
all(is.finite(cont_year$life_w)),
all(is.finite(world_year$world_w))
)
main_lines <- bind_rows(
cont_year %>%
transmute(series = as.character(continent),
year, value = life_w),
world_year %>%
transmute(series = "World", year, value = world_w)
) %>%
mutate(series = factor(series, levels = series_order)) %>%
arrange(series, year)
gap_data <- cont_year %>%
left_join(world_year %>% select(year, world_w), by = "year") %>%
mutate(gap = life_w - world_w)
stopifnot(all(is.finite(gap_data$gap)))
# 4. 跨零填色:插值交点只用于闭合面积,不画作观测
gap_band <- gap_data %>%
group_by(continent) %>%
group_modify(~ {
x <- .x$year
y <- .x$gap
crossing <- which(y[-length(y)] * y[-1L] < 0)
cross_x <- if (length(crossing)) {
x[crossing] -
y[crossing] * (x[crossing + 1L] - x[crossing]) /
(y[crossing + 1L] - y[crossing])
} else numeric()
xx <- sort(unique(c(x, cross_x)))
yy <- stats::approx(x, y, xout = xx)$y
yy[abs(yy) < 1e-10] <- 0
tibble(year = xx, gap = yy)
}) %>%
ungroup()
gap_heads <- gap_data %>%
distinct(continent) %>%
mutate(label = unname(series_zh[as.character(continent)]))
# 5. 共同样式:气泡图同系字体、轴色;补充浅竖向网格
base_theme <- theme_classic(base_family = plot_font, base_size = 11) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.major.x = element_line(
colour = "#EDF1F3", linewidth = 0.18
),
panel.grid.major.y = element_line(
colour = "#E5ECEF", linewidth = 0.22
),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "#738693", linewidth = 0.30),
axis.ticks = element_line(colour = "#738693", linewidth = 0.28),
axis.text = element_text(colour = "#405362", size = 9),
axis.title = element_text(
colour = "#263B4C", face = "bold", size = 10.5
)
)
# 6. 主图:五条系列共用一个图例,色带不进入图例
plot_main <- ggplot() +
geom_ribbon(
data = cont_year,
aes(year, ymin = q25, ymax = q75,
fill = continent, group = continent),
alpha = 0.10, colour = NA, show.legend = FALSE
) +
geom_line(
data = main_lines,
aes(year, value, colour = series,
linetype = series, group = series),
linewidth = 0.87, lineend = "round"
) +
geom_point(
data = main_lines,
aes(year, value, colour = series, shape = series),
size = 1.85, stroke = 0.60
) +
scale_fill_manual(values = continent_cols, guide = "none") +
scale_colour_manual(
name = NULL, values = series_cols,
breaks = series_order, labels = series_zh
) +
scale_linetype_manual(
name = NULL, values = series_linetypes,
breaks = series_order, labels = series_zh
) +
scale_shape_manual(
name = NULL, values = series_shapes,
breaks = series_order, labels = series_zh
) +
guides(
colour = guide_legend(
nrow = 1,
override.aes = list(linewidth = 1.35,
size = 3.20, alpha = 1)
),
linetype = guide_legend(nrow = 1),
shape = guide_legend(nrow = 1)
) +
scale_x_continuous(
breaks = c(1952, 1962, 1972, 1982, 1992, 2002, 2007),
expand = expansion(mult = 0)
) +
scale_y_continuous(
breaks = c(30, 40, 50, 60, 70, 80),
expand = expansion(mult = 0)
) +
coord_cartesian(xlim = c(1950, 2009), ylim = c(28, 85)) +
labs(x = NULL, y = "出生时预期寿命(岁)") +
base_theme +
theme(
legend.position = "top",
legend.justification = "left",
legend.direction = "horizontal",
legend.text = element_text(colour = "#405362", size = 9),
legend.key.width = grid::unit(1.55, "cm"),
legend.key.height = grid::unit(0.70, "cm"),
legend.spacing.x = grid::unit(0.10, "cm"),
plot.margin = margin(7, 13, 6, 10)
)
# 7. 下图:各洲加权均值减 142 国样本参照;虚线为零差值
plot_gaps <- ggplot() +
geom_ribbon(
data = gap_band,
aes(year,
ymin = pmin(gap, 0), ymax = pmax(gap, 0),
fill = continent, group = continent),
alpha = 0.14, colour = NA
) +
geom_hline(
yintercept = 0, linetype = "22",
colour = "#7E8D97", linewidth = 0.48,
show.legend = FALSE
) +
geom_line(
data = gap_data,
aes(year, gap, colour = continent, group = continent),
linewidth = 0.76, lineend = "round"
) +
geom_point(
data = gap_data,
aes(year, gap, colour = continent),
size = 1.45
) +
geom_text(
data = gap_heads,
aes(x = 1979.5, y = 18.3,
label = label, colour = continent),
inherit.aes = FALSE,
family = plot_font, fontface = "bold",
size = 3.45, show.legend = FALSE
) +
facet_wrap(~ continent, ncol = 2) +
scale_colour_manual(values = continent_cols, guide = "none") +
scale_fill_manual(values = continent_cols, guide = "none") +
scale_x_continuous(
breaks = c(1952, 1977, 2007),
expand = expansion(mult = 0)
) +
scale_y_continuous(
breaks = c(-15, 0, 15),
labels = c("−15", "0", "+15"),
expand = expansion(mult = 0)
) +
coord_cartesian(xlim = c(1950, 2009), ylim = c(-20, 20)) +
labs(x = "年份", y = "相对 142 国样本的差值(岁)") +
base_theme +
theme(
legend.position = "none",
strip.text = element_blank(),
strip.background = element_blank(),
panel.spacing = grid::unit(1.5, "lines"),
plot.margin = margin(4, 18, 9, 10)
)
plot_line_advanced <- gridExtra::arrangeGrob(
plot_main, plot_gaps, ncol = 1,
heights = c(1.28, 1)
)
grid::grid.newpage()
grid::grid.draw(plot_line_advanced)排名变化图适合比较同一批对象在不同时点的相对名次。横轴表示时间,纵轴表示排名;每条曲线连接同一对象的各次排名,第 1 名位于顶部。
本例使用 R 自带的 ChickWeight 数据,选取饲料 1 中编号 7
至 12 的六只鸡雏,比较第 4、8、12、16、20
天的体重排名。每个测量日都在这六只鸡雏中按体重由高到低重新排序。
饲料 1 中所选六只鸡雏在五个测量日的体重排名;圆点表示实测记录,起止标签标出体重(克)。数据:R datasets::ChickWeight。
本例从开发者 GitHub 仓库安装固定版本的
ggbump。以下安装命令在 Console
中运行一次。
install.packages(c("ggplot2", "dplyr", "ragg", "systemfonts", "remotes"))
remotes::install_github(
"davidsjoberg/ggbump@fe6d5c7a974d4ae421b295cf81974394df481111",
upgrade = "never"
)作图代码先按日龄计算排名,再通过 group
连接同一鸡雏的观测,并用 scale_y_reverse() 将第 1
名放在上方。本例所选测量日没有同重;若其他数据出现并列,需先确定并列名次的处理方式。
# ChickWeight 排名变化图:同一批鸡雏在五个测量日的组内体重排名。
# 运行后生成 plot_bump;print(plot_bump) 可显示图片。
# 包的安装只需在 Console 运行一次,不写入绘图流程。
# 1. 调用包。ggbump 提供 geom_bump(),用于平滑连接排名位置。
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(ggbump)
})
# 2. 选取同一批对象和相同的测量日。
# Chick 是有序因子,需先转字符再转整数,才能保留原来的鸡雏编号。
# 本例的六只鸡雏均属于饲料 1;排名仅在这六只之间计算。
chick_ids <- 7:12
measurement_days <- c(4, 8, 12, 16, 20)
rank_data <- datasets::ChickWeight %>%
transmute(
chick_id = as.integer(as.character(Chick)),
age_day = Time,
weight_g = weight,
diet = as.character(Diet)
) %>%
filter(chick_id %in% chick_ids, age_day %in% measurement_days)
# 检查每只鸡雏在五个时点均有记录,且没有重复记录。
stopifnot(
nrow(rank_data) == length(chick_ids) * length(measurement_days),
!anyNA(rank_data),
!anyDuplicated(rank_data[c("chick_id", "age_day")]),
all(table(rank_data$chick_id) == length(measurement_days)),
all(table(rank_data$age_day) == length(chick_ids)),
all(rank_data$diet == "1")
)
# 3. 在每个测量日分别排名:体重越大,名次越靠前。
# min_rank() 遇到同重时保留并列,不按编号人为拆分名次。
# 当前数据在各测量日没有同重情况,得到 1—6 名。
rank_data <- rank_data %>%
group_by(age_day) %>%
mutate(rank = min_rank(desc(weight_g))) %>%
ungroup() %>%
arrange(chick_id, age_day) %>%
mutate(chick_key = factor(chick_id, levels = chick_ids))
# 仅在起点和终点标注编号与原始体重,避免图内堆叠文字。
start_labels <- rank_data %>%
filter(age_day == min(measurement_days)) %>%
mutate(label = paste0("鸡雏 ", chick_id, " · ", weight_g, " g"))
end_labels <- rank_data %>%
filter(age_day == max(measurement_days)) %>%
mutate(label = paste0("鸡雏 ", chick_id, " · ", weight_g, " g"))
# 4. 固定每只鸡雏的颜色。取 Okabe–Ito 色盘中的六色,省去亮黄色与黑色。
# 颜色表示对象身份;所有对象使用相同的线宽和点形。
chick_colors <- c(
"7" = "#0072B2", "8" = "#E69F00", "9" = "#009E73",
"10" = "#D55E00", "11" = "#56B4E9", "12" = "#CC79A7"
)
font_family <- if (.Platform$OS.type == "windows") "Microsoft YaHei" else "sans"
# Windows 图形设备有时需先注册中文字体;此设置只影响当前 R 会话。
if (.Platform$OS.type == "windows") {
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(font_family)), font_family))
}
ink <- "#29323A"
# 5. 圆点表示真实测量日的排名;曲线只用于连接这些点。
# smooth 控制连接线形状,不生成新的实测数据,也不估计超越发生的日期。
# scale_y_reverse() 把第 1 名放在顶部;横轴仍使用真实日龄。
plot_bump <- ggplot(
rank_data, aes(x = age_day, y = rank, colour = chick_key, group = chick_key)
) +
geom_hline(
yintercept = 1:6, colour = "#E7EBEE", linewidth = 0.3
) +
geom_bump(smooth = 8, linewidth = 1.25, alpha = 0.88) +
geom_point(aes(fill = chick_key), shape = 21,
size = 3.6, stroke = 0.6, colour = "white") +
geom_text(
data = start_labels, aes(x = 3.35, label = label),
hjust = 1, colour = ink, family = font_family, size = 3.7
) +
geom_text(
data = end_labels, aes(x = 20.65, label = label),
hjust = 0, colour = ink, family = font_family, size = 3.7
) +
scale_colour_manual(values = chick_colors, guide = "none") +
scale_fill_manual(values = chick_colors, guide = "none") +
scale_x_continuous(
breaks = measurement_days,
limits = c(0.2, 24),
expand = expansion(mult = 0)
) +
scale_y_reverse(
breaks = 1:6, limits = c(6.35, 0.65),
expand = expansion(mult = 0)
) +
labs(x = "日龄(天)", y = "组内体重排名") +
theme_classic(base_size = 12, base_family = font_family) +
theme(
axis.text = element_text(colour = ink, size = 11),
axis.title = element_text(colour = ink, size = 12),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
axis.line = element_line(colour = "#77828B", linewidth = 0.3),
axis.ticks = element_line(colour = "#77828B", linewidth = 0.3),
plot.margin = margin(12, 16, 12, 12),
plot.background = element_rect(fill = "white", colour = NA)
)
# 运行 print(plot_bump) 可显示图形;需要保存时可另用 ggsave()。
# 本脚本仅生成对象,供教材的代码块调用,不自动写入文件。
print(plot_bump)动画实现可参见 3.11 动画可视化。资料:R 官方数据说明 · ggbump 项目文档 · AI与R语言实践原文。
1 · 图形与适用问题
柱状图(Bar Chart) 是用于比较分类变量在某一数值指标上的差异的标准图形。它以柱子的长度作为主要编码通道,将数值大小映射为可直接比较的视觉差异,适用于跨类别对比、排序展示与结构呈现。
核心用途
绘图函数
核心 ggplot geom:
geom_col()/geom_bar()
geom_col() 与 geom_bar()
的区别
在 ggplot2 中,柱状图常用两类几何对象,选择取决于数值是否已在数据表中:
geom_col():推荐
当数据表中已经存在要展示的数值列(如
GDP、人口、寿命)时使用:
柱高 = y 变量的实际数值
geom_bar():
当数据为原始明细记录,需要对类别进行计数(频数)时使用:
柱高 = 该类别出现次数
2 · 数据处理
本节使用 gapminder 数据,聚焦 2007
年美洲(Americas)国家,筛选人均 GDP(gdpPercap)最高的 Top
10。
year == 2007 且
continent == "Americas",并按 gdpPercap 取 Top
10(不含并列)。关键技巧:类别变量若不重排,通常会按名称字典序显示。为得到“排名视图”,需要按
gdpPercap 对 country 进行重排(reorder() 或
forcats::fct_reorder())。
代码逻辑:
Step 1 筛选范围(年份 × 大洲)
限定 year == 2007
且 continent == "Americas",得到目标比较集合。
Step 2 提取 Top 10(按人均 GDP 排名)
使用
slice_max(order_by = gdpPercap, n = 10, with_ties = FALSE)
获取人均 GDP 最高的 10 个国家,保证排名集合唯一且可复现。
Step 3 重排类别(避免字典序)
在映射中用
reorder(country, gdpPercap)
将国家按数值从低到高排列;配合水平柱状图即可呈现“从高到低”的视觉排序(y 轴自下而上)。
Step 4 绘制柱状图
用
geom_col() 直接把数据表中的人均 GDP 映射为柱长。
library(tidyverse)
library(gapminder)
data_bar <- gapminder %>%
filter(year == 2007, continent == "Americas") %>%
slice_max(gdpPercap, n = 10, with_ties = FALSE)3 · 最小代码与基础图
数据已筛好,这里只需要设定数值、类别和
geom_col()。国家名按数值重排,最大的柱子显示在最上方。
2007 年美洲样本国家中人均 GDP 排名前十的国家。柱长从零起算,用于比较国家间的数值大小。
读图方法
从图中可直接读出的主要信息
4 · 进阶成图与美化思路
从普通柱状图扩展为堆叠柱状图
进阶 · 美化思路
1. 图形定义
标准柱状图强调“类别—数值”的单一对比。堆叠柱状图(Stacked Bar Chart)在同一根柱子内叠加多个子类别,使图形在一个坐标系中同时表达总量与结构两类信息:
使用要点:堆叠柱状图适合“看总量 + 看组成”;但若目标是精确比较不同主类别间的子类别大小,需要谨慎解读(堆叠基线不同会影响比较)。
2. 场景与目标
本案例希望在一张图中同时回答两类问题:
- 各大洲的样本国家 GDP 估算值合计如何比较(总量比较)?
- 每个大洲的样本合计主要由哪些国家贡献(结构识别)?
3. 数据策略
以 2007 年为截面,先用 pop × gdpPercap 计算各样本国家的
GDP 估算值。这一乘积用于课堂演示,并非官方国民账户
GDP;洲内合计也仅限 gapminder 收录的国家和地区。
若把每个大洲的所有国家都分成独立色块,图例会过长。因此将洲内贡献分为第 1 位国家、第 2 位国家、其余样本国家合计三层。同一大洲使用同一色系,并以深浅区分三层;右侧图例列出前两位国家和其余样本国家类别,图下明细表列出取整后的洲内占比。
进阶 · 成图展示
下面先看进阶成图:上图的柱高是各洲所收录国家的 GDP 估算值之和;下图将各洲自身合计设为 100%,比较洲内组成。右侧图例列出每洲的具体国家与“某洲其他”,柱内百分比只标在占比至少 8% 的色块中。大洋洲在这份数据中只有两国,因此没有“其他”色块。
进阶代码分为国家层估算、洲内排名、三层汇总、绘图与图例几个模块。读图时先比较上图的总量,再用下图、右侧图例和明细表核对组成。
2007 年各洲所收录国家 GDP 估算值的总量与洲内组成。上图柱高为各洲样本合计(万亿美元),下图显示以各洲自身合计为 100% 的组成;右侧图例列出贡献最高的两国及该洲其余样本国家(如有)。估算值为人口乘以人均 GDP,并非官方洲际 GDP。
| 大洲 | 贡献第 1 | 贡献第 2 | 其余样本国家 |
|---|---|---|---|
| Asia | 中国(32%) | 日本(19%) | 亚洲其他(49%) |
| Americas | 美国(67%) | 巴西(9%) | 美洲其他(25%) |
| Europe | 德国(18%) | 英国(14%) | 欧洲其他(68%) |
| Africa | 埃及(19%) | 南非(17%) | 非洲其他(64%) |
| Oceania | 澳大利亚(87%) | 新西兰(13%) | — |
美化思路:每洲保留独立色系,深色突出贡献最大的国家,中间色表示第二位,浅色表示其余样本国家;右侧图例写出实际名称。柱顶标总量,空间足够的色块标洲内占比。先比较上图的柱高,再看下图的内部结构;不同柱中不共用基线的色块不适合精确比较。
# 进阶堆叠柱状图:各洲样本 GDP 估算总量与洲内组成
# 数据口径:Gapminder 2007;GDP 估算值 = 人口 × 人均 GDP
# RStudio 绘图窗较窄时可点击 Zoom;适合约 10.6 × 8.2 英寸的画布
library(dplyr)
library(ggplot2)
library(gapminder)
library(scales)
library(gridExtra)
# 1. 分组顺序、名称与旧版 D3 Category10 洲别基色
year_sel <- 2007L
continent_order <- c("Asia", "Americas", "Europe", "Africa", "Oceania")
continent_zh <- c(
Asia = "亚洲", Americas = "美洲", Europe = "欧洲",
Africa = "非洲", Oceania = "大洋洲"
)
country_zh <- c(
China = "中国", Japan = "日本", `United States` = "美国",
Brazil = "巴西", Germany = "德国", `United Kingdom` = "英国",
Egypt = "埃及", `South Africa` = "南非",
Australia = "澳大利亚", `New Zealand` = "新西兰"
)
base_cols <- c(
Asia = "#1F77B4", Americas = "#FF7F0E",
Europe = "#2CA02C", Africa = "#D62728",
Oceania = "#9467BD"
)
# 同洲三阶:首位最深、次位中等、其余最浅;五洲保留各自色相
mix_col <- function(from, to, weight) {
grDevices::rgb(grDevices::colorRamp(c(from, to))(weight) / 255)
}
shade_triplet <- function(base) {
c(
top1 = mix_col(base, "#182630", 0.35),
top2 = mix_col(base, "#FFFFFF", 0.46),
other = mix_col(base, "#FFFFFF", 0.79)
)
}
segment_cols <- c()
for (ct in continent_order) {
shades <- shade_triplet(base_cols[[ct]])
names(shades) <- paste(ct, names(shades), sep = "__")
segment_cols <- c(segment_cols, shades)
}
# 2. 字体:沿用前一张进阶图的跨系统选择,不下载在线字体
os_name <- Sys.info()[["sysname"]]
font_candidates <- if (identical(os_name, "Windows")) {
c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei")
} else if (identical(os_name, "Darwin")) {
c("PingFang SC", "Heiti SC", "Hiragino Sans GB")
} else {
c("Noto Sans CJK SC", "WenQuanYi Micro Hei")
}
available_fonts <- if (requireNamespace("systemfonts", quietly = TRUE)) {
unique(systemfonts::system_fonts()[["family"]])
} else character()
matched_fonts <- font_candidates[font_candidates %in% available_fonts]
plot_font <- if (length(matched_fonts)) matched_fonts[[1]] else "sans"
if (.Platform$OS.type == "windows" && plot_font != "sans") {
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(plot_font)), plot_font))
}
display_country <- function(x) {
label <- unname(country_zh[x])
label[is.na(label)] <- x[is.na(label)]
label
}
# 3. 每洲按 GDP 估算总量排序;平局按国家名称固定顺序
bar_ranked <- gapminder::gapminder %>%
filter(year == year_sel) %>%
mutate(
continent = factor(as.character(continent), levels = continent_order),
country = as.character(country),
gdp_estimate = pop * gdpPercap
) %>%
group_by(continent) %>%
arrange(desc(gdp_estimate), country, .by_group = TRUE) %>%
mutate(rank_in_continent = row_number()) %>%
ungroup()
stopifnot(
nrow(bar_ranked) == 142L,
all(is.finite(bar_ranked$gdp_estimate)),
all(bar_ranked$gdp_estimate > 0)
)
# 4. 首位、次位、其余国家:上图用绝对值,下图用洲内百分比
bar_parts <- bar_ranked %>%
mutate(
role = case_when(
rank_in_continent == 1L ~ "top1",
rank_in_continent == 2L ~ "top2",
TRUE ~ "other"
),
role = factor(role, levels = c("other", "top2", "top1"))
) %>%
group_by(continent, role) %>%
summarise(
gdp_estimate = sum(gdp_estimate),
n_country = n(),
item = if (as.character(first(role)) == "other") {
"其余样本国家"
} else first(country),
.groups = "drop"
) %>%
group_by(continent) %>%
mutate(
continent_total = sum(gdp_estimate),
share = gdp_estimate / continent_total,
gdp_trillion = gdp_estimate / 1e12,
segment_id = paste(as.character(continent), as.character(role), sep = "__"),
item_zh = if_else(as.character(role) == "other",
paste0(unname(continent_zh[as.character(continent)]), "其他"),
display_country(item)),
label_colour = if_else(as.character(role) == "top1",
"#FFFFFF", "#263B4C")
) %>%
ungroup()
# 固定叠放次序:底部其余、中间次位、顶部首位
segment_levels <- c(
paste(continent_order, "other", sep = "__"),
paste(continent_order, "top2", sep = "__"),
paste(continent_order, "top1", sep = "__")
)
bar_parts <- bar_parts %>%
mutate(segment_id = factor(segment_id, levels = segment_levels))
bar_totals <- bar_parts %>%
distinct(continent, continent_total) %>%
mutate(total_trillion = continent_total / 1e12)
stopifnot(
nrow(bar_totals) == 5L,
all(abs(bar_parts %>% group_by(continent) %>%
summarise(s = sum(share), .groups = "drop") %>% pull(s) - 1) < 1e-10),
!any(bar_parts$continent == "Oceania" & bar_parts$role == "other")
)
# 5. 图例数据:按洲列出实际国家;大洋洲只有两项
legend_data <- bar_parts %>%
mutate(
role_chr = as.character(role),
role_order = match(role_chr, c("top1", "top2", "other")),
legend_group_y = c(14, 11.6, 9.2, 6.8, 4.4)[
match(as.character(continent), continent_order)
],
legend_y = legend_group_y - 0.55 * role_order,
label = item_zh
) %>%
arrange(continent, role_order)
legend_continents <- data.frame(
continent = continent_order,
label = unname(continent_zh[continent_order]),
y = c(14, 11.6, 9.2, 6.8, 4.4),
colour = unname(segment_cols[paste(continent_order, "top1", sep = "__")])
)
# 可在 Console 查看未舍入的洲内占比和国家数量
bar_breakdown <- legend_data %>%
select(continent, role = role_chr, item = item_zh,
n_country, gdp_estimate, share)
# 6. 共同外观:白底、淡水平参考线、细轴线
bar_theme <- theme_classic(base_family = plot_font, base_size = 10) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.major.y = element_line(colour = "#E9EEF0", linewidth = 0.20),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "#788B96", linewidth = 0.28),
axis.ticks = element_line(colour = "#788B96", linewidth = 0.26),
axis.text = element_text(colour = "#405362", size = 9),
axis.title = element_text(colour = "#263B4C", face = "bold", size = 10),
plot.margin = margin(5, 12, 5, 8)
)
fill_scale <- scale_fill_manual(
values = segment_cols,
guide = "none"
)
bar_width <- 0.50 # 两幅图共用;如需更窄,可改为 0.46
# 7A. 绝对总量:同洲色系堆叠;柱顶只标洲合计
plot_total <- ggplot(bar_parts,
aes(x = continent, y = gdp_trillion, fill = segment_id)) +
geom_col(width = bar_width, colour = "white", linewidth = 0.24,
position = position_stack(reverse = TRUE)) +
geom_text(data = bar_totals,
aes(x = continent, y = total_trillion,
label = sprintf("%.1f", total_trillion)),
inherit.aes = FALSE, vjust = -0.50,
colour = "#263B4C", family = plot_font,
fontface = "bold", size = 3.35) +
fill_scale +
scale_x_discrete(labels = continent_zh) +
scale_y_continuous(breaks = c(0, 5, 10, 15, 20),
expand = expansion(mult = c(0, 0.12))) +
labs(x = NULL, y = "样本 GDP 估算合计(万亿美元)") +
bar_theme +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
# 7B. 洲内 100%:同一批数值换成占比;8% 以上的段直接标百分比
plot_share <- ggplot(bar_parts,
aes(x = continent, y = share, fill = segment_id)) +
geom_col(width = bar_width, colour = "white", linewidth = 0.24,
position = position_stack(reverse = TRUE)) +
geom_text(aes(label = if_else(share >= 0.08,
percent(share, accuracy = 1), ""),
colour = label_colour),
position = position_stack(vjust = 0.5, reverse = TRUE),
family = plot_font, fontface = "bold", size = 3.1,
show.legend = FALSE) +
fill_scale +
scale_colour_identity() +
scale_x_discrete(labels = continent_zh) +
scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), labels = label_percent(),
expand = expansion(mult = c(0, 0))) +
labs(x = NULL, y = "洲内占比") +
bar_theme +
theme(plot.margin = margin(4, 12, 5, 8))
# 8. 右侧分组图例:洲名下只列实际国家和该洲其他
bar_legend <- ggplot(legend_data) +
geom_text(data = legend_continents,
aes(x = 0.08, y = y, label = label, colour = colour),
inherit.aes = FALSE, hjust = 0,
family = plot_font, fontface = "bold", size = 3.55) +
geom_rect(aes(xmin = 0.08, xmax = 0.20,
ymin = legend_y - 0.16, ymax = legend_y + 0.16,
fill = segment_id),
colour = "#A5B2BA", linewidth = 0.17) +
geom_text(aes(x = 0.28, y = legend_y, label = label),
hjust = 0, family = plot_font, size = 3.2,
colour = "#405362") +
scale_fill_manual(values = segment_cols, guide = "none") +
scale_colour_identity() +
scale_x_continuous(limits = c(0, 1), expand = expansion(mult = 0)) +
scale_y_continuous(limits = c(0, 17.4), expand = expansion(mult = 0)) +
theme_void(base_family = plot_font) +
theme(plot.margin = margin(4, 5, 4, 12))
# 9. 对齐上下两图的绘图区,图例置于右侧;RStudio 中直接显示
grob_total <- ggplotGrob(plot_total)
grob_share <- ggplotGrob(plot_share)
shared_widths <- grid::unit.pmax(grob_total$widths, grob_share$widths)
grob_total$widths <- shared_widths
grob_share$widths <- shared_widths
bar_main <- gridExtra::arrangeGrob(
grob_total, grob_share,
ncol = 1, heights = c(1.42, 1)
)
plot_bar_advanced <- gridExtra::arrangeGrob(
bar_main, bar_legend,
ncol = 2, widths = c(3.65, 1)
)
grid::grid.newpage()
grid::grid.draw(plot_bar_advanced)
# 10. 图下明细表沿用同一批汇总值,给出各洲具体名称和占比
continent_breakdown <- bar_parts %>%
mutate(
cell = paste0(item_zh, "(", scales::percent(share, accuracy = 1), ")"),
role_chr = factor(as.character(role), levels = c("top1", "top2", "other"))
) %>%
select(continent, role_chr, cell) %>%
tidyr::pivot_wider(
names_from = role_chr,
values_from = cell,
values_fill = "—"
) %>%
select(continent, top1, top2, other) %>%
arrange(continent)1 · 图形与适用问题
直方图(Histogram) 用于刻画单个连续变量的分布形态。它将连续取值按给定组距划分为若干区间(Bins),并统计每个区间内的观测频数,从而以“柱形轮廓”近似呈现数据的整体结构与集中趋势。
核心用途
与柱状图的区别
geom_histogram()):geom_col()
/ geom_bar()):核心函数
在 ggplot2 中,直方图对应的几何对象函数为:
geom_histogram()
以下用 2007 年各国预期寿命示范。先准备单年数据,再出基础图和进阶图。
2 · 数据处理
先筛选 2007 年,每行代表一个国家或地区;下面两张图均沿用这个截面。
3 · 最小代码与基础图
横轴放连续变量,geom_histogram()
统计区间内的观测数;binwidth = 2 表示每组宽 2 年。
ggplot(data_hist, aes(x = lifeExp)) +
geom_histogram(binwidth = 2) +
labs(x = "预期寿命(岁)", y = "国家或地区数")2007 年各国预期寿命的频数分布;每个国家或地区计一条观测,组距为 2 年。
基础图怎么读? 横轴是预期寿命区间,纵轴是该区间的国家或地区数;柱越高,该区间的观测越多。改变组距会改变柱形轮廓,细小峰谷应结合组距检查。这里每个国家或地区计一次,不按人口加权。
4 · 进阶成图与美化思路
进阶图从 2007 年数据中选取亚洲、美洲、欧洲和非洲,共 140 个国家或地区,以 2×2 分面比较预期寿命分布。每个分面叠加密度直方图与 KDE 曲线,并标出均值、中位数和第 1、3 四分位数。四洲共用 35—90 岁的横轴与相同的密度纵轴,便于比较分布形态;密度柱高不能用来比较各洲国家数量。大洋洲只有 2 条记录,未纳入此处分布比较。
图上的三个检查点:分面标题标明各洲样本量;右侧图例以虚线、实线和点线分别说明均值、中位数与 Q1/Q3,底部深色横棒连接 Q1—Q3,白心点标出中位数;每道底部细刻度对应一个国家或地区。直方图组距为 4 年,KDE 带宽为 3 年。这两项设定都会影响轮廓细节,读图时应结合原始刻度判断。
2007 年亚洲、美洲、欧洲和非洲国家或地区的预期寿命分布。四个分面叠加 4 年组距的密度直方图与带宽 3 年的 KDE 曲线;右侧图例区分均值、中位数及 Q1/Q3 位置线,底部短刻度为原始观测,深色横棒和白心点分别表示四分位区间和中位数。各国家或地区等权;大洋洲仅 2 条记录,未纳入分布比较。
suppressPackageStartupMessages({
library(tidyverse)
library(gapminder)
})
# 1. 数据:2007 年每个国家或地区计一次,不按人口加权。
# 大洋洲仅有 2 条记录,暂不纳入分布比较。
continent_order <- c("Asia", "Americas", "Europe", "Africa")
continent_zh <- c(Asia = "亚洲", Americas = "美洲",
Europe = "欧洲", Africa = "非洲")
fill_cols <- c(Asia = "#59977C", Americas = "#5D88AD",
Europe = "#977AA8", Africa = "#B76B72")
line_cols <- c(Asia = "#397D63", Americas = "#376F9A",
Europe = "#795A94", Africa = "#A9515E")
hist_data <- gapminder %>%
filter(year == 2007, as.character(continent) %in% continent_order) %>%
mutate(continent = factor(as.character(continent), levels = continent_order))
stopifnot(nrow(hist_data) == 140L)
# 2. 每洲描述性统计:每个国家或地区的权重相同。
hist_stats <- hist_data %>%
group_by(continent) %>%
summarise(
n = n(),
mean_life = mean(lifeExp),
median_life = median(lifeExp),
q1 = unname(quantile(lifeExp, 0.25)),
q3 = unname(quantile(lifeExp, 0.75)),
.groups = "drop"
) %>%
mutate(
stat_label = sprintf(
"均值 %.1f · 中位数 %.1f · Q1 %.1f · Q3 %.1f",
mean_life, median_life, q1, q3
)
)
facet_labels <- setNames(
paste0(continent_zh[as.character(hist_stats$continent)],
" · n = ", hist_stats$n),
as.character(hist_stats$continent)
)
# 3. KDE:各洲采用相同的 3 年带宽。
# 直方图纵轴为密度,不能以柱高比较各洲国家数量。
density_lines <- hist_data %>%
group_by(continent) %>%
group_modify(~ {
d <- stats::density(.x$lifeExp, bw = 3, from = 35, to = 90, n = 512)
tibble(lifeExp = d$x, density = d$y)
}) %>%
ungroup()
# 4. 统计位置线:Q1 和 Q3 共用一项图例,但在图内各画一条线。
stat_vlines <- bind_rows(
hist_stats %>%
transmute(continent, stat_type = "均值",
stat_x = mean_life, y_end = 0.118),
hist_stats %>%
transmute(continent, stat_type = "中位数",
stat_x = median_life, y_end = 0.118),
hist_stats %>%
transmute(continent, stat_type = "Q1 / Q3",
stat_x = q1, y_end = 0.052),
hist_stats %>%
transmute(continent, stat_type = "Q1 / Q3",
stat_x = q3, y_end = 0.052)
) %>%
mutate(
stat_type = factor(stat_type,
levels = c("均值", "中位数", "Q1 / Q3"))
)
# 5. 风格:统一字体、较深线条和克制的辅助网格。
os_name <- Sys.info()[["sysname"]]
font_candidates <- if (identical(os_name, "Windows")) {
c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei")
} else if (identical(os_name, "Darwin")) {
c("PingFang SC", "Heiti SC", "Hiragino Sans GB")
} else {
c("Noto Sans CJK SC", "WenQuanYi Micro Hei")
}
available_fonts <- if (requireNamespace("systemfonts", quietly = TRUE)) {
unique(systemfonts::system_fonts()[["family"]])
} else character()
matched_fonts <- font_candidates[font_candidates %in% available_fonts]
plot_font <- if (length(matched_fonts)) matched_fonts[[1]] else "sans"
if (.Platform$OS.type == "windows" && plot_font != "sans") {
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(plot_font)), plot_font))
}
ink <- "#293E4A"
muted <- "#61717C"
plot_hist_advanced <- ggplot(hist_data, aes(x = lifeExp)) +
# 4 年组距的密度直方图。
geom_histogram(
aes(y = after_stat(density), fill = continent),
binwidth = 4, boundary = 0, closed = "left",
colour = "white", linewidth = 0.25, alpha = 0.68
) +
# 同一带宽的 KDE 曲线;颜色与柱形保持一致。
geom_line(
data = density_lines,
aes(x = lifeExp, y = density, colour = continent),
inherit.aes = FALSE,
linewidth = 1.12, lineend = "round",
show.legend = FALSE
) +
# 均值、中位数、Q1/Q3 的位置;线型自动形成右侧图例。
geom_segment(
data = stat_vlines,
aes(x = stat_x, xend = stat_x,
y = 0, yend = y_end, linetype = stat_type),
inherit.aes = FALSE,
colour = ink, linewidth = 0.64, alpha = 0.86
) +
geom_hline(yintercept = 0, colour = "#8D9BA2", linewidth = 0.34) +
# 每条底部短刻度表示一个国家或地区。
geom_point(
aes(y = -0.006, colour = continent),
shape = "|", size = 2.45, alpha = 0.72,
show.legend = FALSE
) +
# Q1—Q3 横棒,白心点为中位数。
geom_segment(
data = hist_stats,
aes(x = q1, xend = q3, y = -0.019, yend = -0.019),
inherit.aes = FALSE,
colour = ink, linewidth = 1.75, lineend = "round"
) +
geom_point(
data = hist_stats,
aes(x = median_life, y = -0.019),
inherit.aes = FALSE,
shape = 21, size = 2.5,
fill = "white", colour = ink, stroke = 0.68
) +
# 将数值放在固定位置,避免邻近统计量的标签互相遮挡。
geom_text(
data = hist_stats,
aes(x = 36.5, y = -0.044, label = stat_label),
inherit.aes = FALSE,
hjust = 0, vjust = 0.5,
family = plot_font, size = 2.85,
colour = "#344B59"
) +
facet_wrap(
~ continent, ncol = 2, scales = "fixed",
labeller = as_labeller(facet_labels)
) +
scale_fill_manual(values = fill_cols, guide = "none") +
scale_colour_manual(values = line_cols, guide = "none") +
scale_linetype_manual(
name = "统计位置",
values = c("均值" = "dashed",
"中位数" = "solid",
"Q1 / Q3" = "dotted"),
breaks = c("均值", "中位数", "Q1 / Q3")
) +
scale_x_continuous(
breaks = seq(40, 90, 10),
expand = c(0, 0)
) +
scale_y_continuous(
breaks = c(0, 0.04, 0.08, 0.12),
labels = c("0", "0.04", "0.08", "0.12"),
expand = c(0, 0)
) +
coord_cartesian(
xlim = c(35, 90),
ylim = c(-0.057, 0.125)
) +
labs(
x = "出生时预期寿命(岁)",
y = "密度"
) +
guides(
linetype = guide_legend(
override.aes = list(linewidth = 1.2, alpha = 1)
)
) +
theme_classic(base_family = plot_font, base_size = 10) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.major.y = element_line(colour = "#E7EDF0", linewidth = 0.25),
panel.grid.major.x = element_line(colour = "#F0F3F5", linewidth = 0.18),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "#71818A", linewidth = 0.30),
axis.ticks = element_line(colour = "#71818A", linewidth = 0.27),
axis.text = element_text(colour = "#435360", size = 8.5),
axis.title = element_text(colour = ink, face = "bold", size = 10),
strip.background = element_rect(fill = "#EDF2F4", colour = NA),
strip.text = element_text(colour = ink, face = "bold", size = 9.5),
panel.spacing = grid::unit(0.68, "lines"),
legend.position = "right",
legend.title = element_text(colour = ink, face = "bold", size = 9.5),
legend.text = element_text(colour = "#344B59", face = "bold", size = 9),
legend.key.width = grid::unit(1.1, "cm"),
legend.key.height = grid::unit(0.60, "cm"),
legend.spacing.y = grid::unit(0.14, "cm"),
plot.margin = margin(8, 11, 8, 8)
)
print(plot_hist_advanced)KDE 如何得到一条曲线? 在每个观测值的位置放上一条小的核函数曲线,把它们相加并按样本量和带宽调整,就得到图中的总体密度曲线。调大带宽会让曲线更平滑;调小带宽会保留更多局部起伏,但小峰不一定对应真实的人群分组。
上面的互动组件使用 0—100
的示意数值来展示计算过程,和本节 gapminder
的预期寿命直方图不是同一份数据。
1 · 图形与适用问题
雷达图(Radar Chart)把多个指标分别放在从中心向外的轴上,用折线连接一个对象在各轴上的位置。它适合指标数量不多、指标含义清楚、需要快速比较结构轮廓的情境,例如比较数辆车在油耗、马力、车重等指标上的相对位置。各对象必须使用同一组指标与同一套尺度。
本例使用 mtcars 中的五个变量:mpg
为每加仑行驶英里数,hp 为马力,qsec
为行驶四分之一英里所需秒数,wt
为车重(千磅),drat
为后桥传动比。这些变量单位不同,先分别缩放到
0—1。越靠外只表示原始数值越大;例如 qsec
越大反而表示完成相同距离所用时间越长,不能把整个多边形面积当作综合评分。改变指标顺序也会改变多边形形状。
基础图使用
fmsb::radarchart(),便于先用一个函数认识雷达图;进阶图再用
ggplot2
计算轴上的二维坐标,并通过分面比较多款车型。fmsb 属于 R
的基础绘图体系,不沿用 ggplot() + geom_*()
语法。首次使用可在 Console 运行
install.packages("fmsb")。
2 · 数据处理
先按 mtcars 全部 32 款车型的每个指标分别做 min–max
缩放,再取 Mazda RX4。这样得到的 0—1
是样本内相对位置,并非车型的绝对性能。前两行全 1 和全 0
是 fmsb 固定轴界所需的数据格式。
library(tidyverse)
library(fmsb)
vars <- c("mpg", "hp", "qsec", "wt", "drat")
radar_scores <- mtcars %>%
rownames_to_column("car") %>%
select(car, all_of(vars)) %>%
pivot_longer(-car, names_to = "var", values_to = "raw") %>%
group_by(var) %>%
mutate(score = (raw - min(raw)) / (max(raw) - min(raw))) %>%
ungroup() %>%
filter(car == "Mazda RX4") %>%
select(var, score) %>%
pivot_wider(names_from = var, values_from = score)
radar_in <- as.data.frame(rbind(
rep(1, length(vars)),
rep(0, length(vars)),
unlist(radar_scores[vars], use.names = FALSE)
))
names(radar_in) <- vars3 · 最小代码与基础图
数据整理好后,绘图只需调用一个函数。为让形状可见,这里只加轮廓色、透明填充和轴刻度;这些设置不参与标准化计算。
fmsb::radarchart(
radar_in,
pcol = "#356F8E",
pfcol = scales::alpha("#356F8E", 0.18),
plwd = 2,
axistype = 1
)Mazda RX4 在 mtcars 样本内五项指标的 0—1 相对位置;越靠外仅表示原始数值越大。
基础图怎么读? 沿每条轴看 Mazda RX4 的顶点。0
对应该指标在这 32 款车中的最小值,1
对应最大值;顶点越靠外,表示该变量原始数值越大。比较不同车型时,需要逐轴比较,不能把多边形面积或“凹陷”直接解释为总成绩与短板。尤其要留意
mpg、qsec、wt
的方向含义各不相同。
4 · 进阶成图与美化思路
进阶图进一步比较这 32 款车中 mpg 与 hp
各自最高、最低的车型。四张雷达图分别展示一款车的五项标准化指标,并以灰色虚线给出全样本的逐指标中位数。若出现并列极值,代码按车型名称排序后取第一款。
美化与读图思路。 四个分面共用五条轴与 0—1
半径,四圈五边形网格依次为
0.25、0.50、0.75、1.00。四款车分别以蓝、橙、绿、紫色实线、半透明填充和空心点显示;灰色虚线连接
32
款车各指标的中位数,只用于逐轴参照,不代表一款真实车型。沿同一轴比较车型顶点与参考线:mpg
大通常表示单位燃油行驶更远,qsec
大表示完成四分之一英里所需时间更长,wt
大表示车更重。向外延伸只说明原始数值较大,不能把轮廓面积当成综合成绩。
进阶图继续用 ggplot2
组合网格、参考线和车型轮廓;完整代码收在图下折叠区。
mtcars 中按 mpg 与 hp 极值选出的四款车。各分面蓝、橙、绿、紫轮廓分别表示该车的五项样本内标准化值;灰色虚线为 32 款车逐指标中位数,四圈网格从内到外为 0.25、0.50、0.75、1.00。越靠外仅表示该指标原始值较大。
suppressPackageStartupMessages({
library(tidyverse)
})
# 1. 选车 ----------------------------------------------------
car_tbl <- mtcars %>%
tibble::rownames_to_column("car")
selected_cars <- bind_rows(
car_tbl %>%
arrange(desc(mpg), car) %>%
slice(1) %>%
mutate(choice = "mpg 最高"),
car_tbl %>%
arrange(mpg, car) %>%
slice(1) %>%
mutate(choice = "mpg 最低"),
car_tbl %>%
arrange(desc(hp), car) %>%
slice(1) %>%
mutate(choice = "hp 最高"),
car_tbl %>%
arrange(hp, car) %>%
slice(1) %>%
mutate(choice = "hp 最低")
) %>%
distinct(car, .keep_all = TRUE) %>%
mutate(
panel = paste0(choice, "\n", car)
)
stopifnot(n_distinct(selected_cars$car) == 4L)
panel_order <- selected_cars$panel
# 四款车型分别使用四种颜色
car_cols <- setNames(
c(
"#2E86AB", # 蓝
"#D16D3B", # 橙
"#4C956C", # 绿
"#8E63A9" # 紫
),
selected_cars$car
)
median_colour <- "#7A8D98"
# 2. 指标与坐标 ----------------------------------------------
metric_order <- c("mpg", "hp", "qsec", "wt", "drat")
axis_info <- tibble(
metric = metric_order,
axis_index = seq_along(metric_order),
angle = pi / 2 -
2 * pi * (axis_index - 1) / length(metric_order),
axis_label = c(
"mpg\n里程/加仑",
"hp\n马力",
"qsec\n1/4 英里用时",
"wt\n车重(千磅)",
"drat\n传动比"
)
) %>%
mutate(
xend = cos(angle),
yend = sin(angle),
label_x = 1.18 * xend,
label_y = 1.18 * yend,
label_hjust = case_when(
abs(xend) < 0.10 ~ 0.5,
xend > 0 ~ 0,
TRUE ~ 1
)
)
# 3. 标准化得分 ----------------------------------------------
# 各指标按全部 32 款车的最小值和最大值缩放到 0—1。
# 越靠外表示原始数值越大,不表示综合性能越好。
all_scores <- car_tbl %>%
select(car, all_of(metric_order)) %>%
pivot_longer(
-car,
names_to = "metric",
values_to = "raw"
) %>%
group_by(metric) %>%
mutate(
score = (raw - min(raw)) /
(max(raw) - min(raw))
) %>%
ungroup() %>%
left_join(
axis_info %>%
select(metric, axis_index, angle),
by = "metric"
)
selected_scores <- all_scores %>%
inner_join(
selected_cars %>% select(car, panel, choice),
by = "car"
) %>%
mutate(
panel = factor(panel, levels = panel_order),
x = score * cos(angle),
y = score * sin(angle)
) %>%
arrange(panel, axis_index)
# 显式闭合轮廓
selected_path <- selected_scores %>%
group_by(panel, car, choice) %>%
group_modify(~ {
d <- arrange(.x, axis_index)
bind_rows(d, d[1, ])
}) %>%
ungroup()
# 4. 样本中位数参考轮廓 --------------------------------------
reference <- all_scores %>%
group_by(metric, axis_index, angle) %>%
summarise(
med = median(score),
.groups = "drop"
) %>%
arrange(axis_index) %>%
mutate(
x_med = med * cos(angle),
y_med = med * sin(angle)
)
median_path <- bind_rows(
reference,
reference[1, ]
)
# 5. 网格与半径刻度 ------------------------------------------
grid_rings <- tidyr::crossing(
radius = c(0.25, 0.50, 0.75, 1.00),
axis_info
) %>%
arrange(radius, axis_index) %>%
mutate(
x = radius * cos(angle),
y = radius * sin(angle)
)
ring_labels <- tibble(
x = 0.04,
y = c(0.25, 0.50, 0.75, 1.00),
label = c("0.25", "0.50", "0.75", "1.00")
)
# 6. 中文字体 ------------------------------------------------
font_candidates <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Heiti SC", "Arial Unicode MS"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
font_available <- tryCatch(
unique(systemfonts::system_fonts()$family),
error = function(e) character()
)
font_match <- font_candidates[
font_candidates %in% font_available
]
plot_font <- if (length(font_match)) {
font_match[1]
} else {
"sans"
}
# 7. 作图 ----------------------------------------------------
radar_refined <- ggplot() +
# 五边形网格
geom_polygon(
data = grid_rings,
aes(x, y, group = radius),
fill = NA,
colour = "#D9E3E8",
linewidth = 0.34,
show.legend = FALSE
) +
# 中心轴
geom_segment(
data = axis_info,
aes(
x = 0, y = 0,
xend = xend, yend = yend
),
colour = "#D9E3E8",
linewidth = 0.34,
show.legend = FALSE
) +
# 半径标尺
geom_text(
data = ring_labels,
aes(x, y, label = label),
family = plot_font,
size = 2.55,
colour = "#8FA0AA",
hjust = 0,
show.legend = FALSE
) +
# 样本中位数:
# 将 linetype 放入 aes(),单独生成参考线图例。
geom_path(
data = median_path,
aes(
x = x_med,
y = y_med,
group = 1,
linetype = "样本中位数"
),
colour = median_colour,
linewidth = 0.82,
lineend = "round",
key_glyph = "path",
show.legend = TRUE
) +
# 车型填充面
geom_polygon(
data = selected_path,
aes(
x, y,
group = panel,
fill = car
),
alpha = 0.16,
colour = NA,
show.legend = FALSE
) +
# 车型轮廓
geom_path(
data = selected_path,
aes(
x, y,
group = panel,
colour = car
),
linewidth = 1.12,
lineend = "round",
show.legend = FALSE
) +
# 车型顶点
geom_point(
data = selected_scores,
aes(x, y, colour = car),
shape = 21,
size = 2.65,
stroke = 0.92,
fill = "white",
show.legend = FALSE
) +
# 轴标签
geom_text(
data = axis_info,
aes(
x = label_x,
y = label_y,
label = axis_label,
hjust = label_hjust
),
family = plot_font,
size = 3.25,
lineheight = 0.92,
colour = "#405463",
show.legend = FALSE
) +
facet_wrap(
~ panel,
ncol = 2
) +
coord_equal(
xlim = c(-1.38, 1.38),
ylim = c(-1.22, 1.28),
expand = FALSE,
clip = "off"
) +
# 车型名称已经写在分面标题中,不重复生成颜色图例
scale_colour_manual(
values = car_cols,
guide = "none"
) +
scale_fill_manual(
values = car_cols,
guide = "none"
) +
# 独立的中位数虚线图例
scale_linetype_manual(
name = NULL,
values = c("样本中位数" = "22"),
breaks = "样本中位数"
) +
guides(
linetype = guide_legend(
nrow = 1,
override.aes = list(
colour = median_colour,
linewidth = 0.82,
alpha = 1
)
)
) +
theme_void(
base_family = plot_font,
base_size = 11
) +
theme(
plot.background = element_rect(
fill = "white",
colour = NA
),
panel.background = element_rect(
fill = "white",
colour = NA
),
strip.background = element_rect(
fill = "#F3F6F8",
colour = NA
),
strip.text = element_text(
colour = "#263E4C",
face = "bold",
size = 10.2,
lineheight = 1.04
),
panel.spacing = grid::unit(1.0, "lines"),
# 底部只显示一项中位数图例
legend.position = "bottom",
legend.direction = "horizontal",
legend.title = element_blank(),
legend.text = element_text(
colour = "#3C5260",
size = 10
),
legend.key.width = grid::unit(1.2, "cm"),
legend.key.height = grid::unit(0.4, "cm"),
legend.background = element_rect(
fill = "white",
colour = NA
),
legend.margin = margin(3, 0, 0, 0),
plot.margin = margin(10, 22, 12, 22)
)
radar_refined1 · 图形与适用问题
箱线图(Boxplot)用少量统计位置概括一组连续数值。当手头有一个类别变量和一个连续变量,例如大洲与各国预期寿命,且希望比较各组的中心位置、分散程度和远离主体的观测时,箱线图比逐组列出统计表更便于横向查看。它不显示分布的全部细节;样本少的组尤其需要同时查看原始点。
箱体从第 25 百分位数(\(Q_1\))到第 75 百分位数(\(Q_3\)),表示中间约一半的观测;箱内横线为中位数。四分位距 \(IQR=Q_3-Q_1\)。默认胡须的端点是落在 \([Q_1-1.5IQR,\;Q_3+1.5IQR]\) 范围内的最远实际观测值,并不一定等于这两个计算边界。范围外的观测单独画成点;它们是核查线索,不能仅凭图形判定为数据错误。
2 · 数据处理
基础图与进阶图都使用 2007 年数据,每行代表一个国家或地区。基础图展示五洲的 142 条记录;进阶图为便于比较分布形态,排除仅有两国的大洋洲,展示其余四洲的 140 条记录。
library(tidyverse)
library(gapminder)
data_box <- gapminder %>%
filter(year == 2007)
continent_labels <- c(
Africa = "非洲", Americas = "美洲", Asia = "亚洲",
Europe = "欧洲", Oceania = "大洋洲"
)3 · 最小代码与基础图
continent 是分组,lifeExp
是组内要概括的连续变量;一个 geom_boxplot()
就能得到各大洲的基础箱线图。
ggplot(data_box, aes(x = continent, y = lifeExp)) +
geom_boxplot() +
scale_x_discrete(labels = continent_labels) +
labs(x = "大洲", y = "预期寿命(岁)")2007 年各大洲国家预期寿命的箱线图;各国等权,大洋洲仅有两条观测。
基础图怎么读? 先比较箱内横线的位置,再比较箱体高度与胡须延伸范围。这里的每个观测是一个国家,不是一个居民;图形反映的是各国数值的分布,不代表按人口加权的居民寿命分布。大洋洲每年只有澳大利亚和新西兰两条记录,箱体与胡须不足以展示稳定的分布形态,应查看两个原始值。
4 · 进阶成图与美化思路
进阶图使用非洲、亚洲、美洲和欧洲的 140 个国家,将半小提琴密度、箱线图与国家点阵叠在同一数值轴上。彩色轮廓显示各洲国家数值较集中的区间;箱体概括中间一半的国家,方形点则保留每个国家的观测。大洋洲仅有两国,不纳入这张密度比较图。
美化思路。 四洲使用与前面图表一致的低饱和配色;方形点按预期寿命着色,白色细边帮助辨认重叠观测。红菱形标出各洲国家的均值,白色箱体的中线标出中位数。图顶端仅保留校正后达到显著水平的组间比较括号,避免不显著比较挤占图面。
2007 年四洲国家预期寿命的云雨图与组间比较
以下三段按顺序在同一 R 会话运行;需预先安装
tidyverse、gapminder、ggdist、rstatix、systemfonts
和 ragg。统计结果由代码计算,图上的括号只负责显示结果。
① 数据准备:筛选四洲并设置标签与颜色
library(tidyverse)
library(gapminder)
library(ggdist)
# 四洲共用 2007 年的国家截面,并固定绘图顺序
rain_order <- c("Africa", "Asia", "Americas", "Europe")
rain_cols <- c(
Africa = "#B76B72", Asia = "#59977C",
Americas = "#5D88AD", Europe = "#977AA8"
)
rain_names <- c(
Africa = "非洲", Asia = "亚洲",
Americas = "美洲", Europe = "欧洲"
)
rain_df <- gapminder %>%
filter(year == 2007, continent %in% rain_order) %>%
mutate(
continent = factor(as.character(continent), levels = rev(rain_order)),
group_id = as.numeric(continent)
)
rain_n <- rain_df %>% count(continent)
rain_labels <- setNames(
paste0(rain_names[as.character(rain_n$continent)], "\nn = ", rain_n$n),
as.character(rain_n$continent)
)② 显著性检验:总体比较与 Holm 校正后的两两比较
# 2. 组间检验:Kruskal–Wallis 总体检验与 Dunn 两两比较
# 六组两两比较统一做 Holm 校正;国家为观测单位,不作居民层面或因果推断。
kw_test <- kruskal.test(lifeExp ~ continent, data = rain_df)
dunn_results <- rstatix::dunn_test(
rain_df, lifeExp ~ continent, p.adjust.method = "holm"
) %>%
mutate(
pair_key = paste(pmin(group1, group2), pmax(group1, group2), sep = "__"),
stars = case_when(
p.adj < 0.001 ~ "***",
p.adj < 0.01 ~ "**",
p.adj < 0.05 ~ "*",
TRUE ~ "ns"
)
)
# 只绘制显著比较;亚洲—美洲的结果保留在检验表和图注中。
# 括号位置只控制排版,不参与统计计算。
bracket_layout <- tribble(
~group1, ~group2, ~bar_at,
"Americas", "Europe", 89.0,
"Africa", "Asia", 91.8,
"Asia", "Europe", 94.6,
"Africa", "Americas", 97.4,
"Africa", "Europe", 100.2
) %>%
mutate(pair_key = paste(pmin(group1, group2), pmax(group1, group2), sep = "__"))
brackets <- bracket_layout %>%
left_join(
dunn_results %>% select(pair_key, p.adj, stars),
by = "pair_key"
) %>%
filter(stars != "ns") %>%
mutate(
x1 = match(group1, levels(rain_df$continent)),
x2 = match(group2, levels(rain_df$continent)),
xmid = (x1 + x2) / 2
)
stopifnot(nrow(brackets) == 5, !anyNA(brackets$p.adj))③ 绘图:半小提琴、箱线、点阵与比较括号
# 3. 字体:使用系统已有的中文无衬线字体,不下载字体
font_candidates <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Hiragino Sans GB", "Heiti SC"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
installed_fonts <- tryCatch(
unique(systemfonts::system_fonts()$family),
error = function(e) character()
)
font_match <- font_candidates[font_candidates %in% installed_fonts]
plot_font <- if (length(font_match)) font_match[[1]] else "sans"
# 4. 进阶云雨图:半小提琴 → 箱线 → 均值 → 国家点阵 → 比较括号
raincloud_plot <- ggplot(
rain_df,
aes(x = group_id, y = lifeExp, group = continent)
) +
ggdist::stat_halfeye(
aes(colour = continent, fill = after_scale(scales::alpha(colour, 0.56))),
# 核密度尾部自然收至零;曲线延伸不代表存在额外国家观测。
density = "unbounded",
trim = FALSE,
adjust = 0.8,
width = 0.50,
justification = -0.22,
scale = 0.66,
normalize = "groups",
point_interval = NULL,
.width = NULL,
slab_linewidth = 0.32,
show.legend = FALSE
) +
geom_boxplot(
width = 0.16,
fill = "white",
colour = "#243640",
linewidth = 0.50,
outlier.shape = NA,
position = position_nudge(x = -0.04),
show.legend = FALSE
) +
stat_summary(
fun = mean, geom = "point",
# 23 是可填充的菱形;红色实心与白色箱体和中位数线区分。
shape = 23, size = 3.15,
fill = "#B84855", colour = "#772A36", stroke = 0.65,
position = position_nudge(x = -0.04),
show.legend = FALSE
) +
ggdist::geom_dots(
aes(fill = lifeExp),
# 方形点形成点阵,白色细边将相互接触的国家点隔开。
side = "left", justification = 1.34,
binwidth = 0.85, overflow = "compress",
scale = 0.41, dotsize = 0.90, stackratio = 0.88,
shape = 22, colour = "white", linewidth = 0.18,
alpha = 0.94
) +
geom_segment(
data = brackets,
aes(x = x1, xend = x2, y = bar_at, yend = bar_at),
inherit.aes = FALSE,
linewidth = 0.48, colour = "#50616B",
show.legend = FALSE
) +
geom_segment(
data = brackets,
aes(x = x1, xend = x1, y = bar_at - 0.60, yend = bar_at),
inherit.aes = FALSE,
linewidth = 0.48, colour = "#50616B",
show.legend = FALSE
) +
geom_segment(
data = brackets,
aes(x = x2, xend = x2, y = bar_at - 0.60, yend = bar_at),
inherit.aes = FALSE,
linewidth = 0.48, colour = "#50616B",
show.legend = FALSE
) +
# 星号放在括号外侧,留足与自身及相邻比较线的距离。
geom_text(
data = brackets,
aes(x = xmid, y = bar_at + 1.10, label = stars),
inherit.aes = FALSE,
family = plot_font, fontface = "bold",
colour = "#A33745", size = 3.5,
show.legend = FALSE
) +
scale_colour_manual(values = rain_cols, guide = "none") +
scale_fill_viridis_c(
option = "D", direction = -1, begin = 0.07, end = 0.90,
name = "预期寿命(岁)",
breaks = c(40, 55, 70, 80),
guide = guide_colourbar(
direction = "horizontal",
title.position = "top",
label.position = "bottom",
barwidth = grid::unit(38, "mm"),
barheight = grid::unit(3.4, "mm")
)
) +
scale_x_continuous(
breaks = seq_along(levels(rain_df$continent)),
labels = unname(rain_labels[levels(rain_df$continent)]),
expand = expansion(add = c(0.26, 0.55))
) +
scale_y_continuous(
breaks = seq(30, 90, 10),
expand = expansion(mult = c(0.01, 0.01))
) +
coord_flip(clip = "off") +
labs(
x = NULL,
y = "出生时预期寿命(岁)",
caption = paste0(
"每点一国(国家等权);红菱形为均值,箱体中线为中位数;大洋洲 2 国未纳入。\n",
"总体 Kruskal–Wallis P = ",
format(kw_test$p.value, scientific = TRUE, digits = 3),
";两两 Dunn 检验经 Holm 校正。\n",
"仅标显著比较:* P < 0.05,** P < 0.01,*** P < 0.001;亚洲—美洲未达显著(P = ",
format(dunn_results$p.adj[dunn_results$pair_key == "Americas__Asia"], digits = 3),
")。"
)
) +
theme_classic(base_family = plot_font, base_size = 11) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.major.y = element_line(colour = "#E5EBEE", linewidth = 0.28,
linetype = "22"),
panel.grid.major.x = element_line(colour = "#F1F4F5", linewidth = 0.20),
panel.grid.minor = element_blank(),
axis.line.y = element_line(colour = "#334650", linewidth = 0.48),
axis.line.x = element_line(colour = "#334650", linewidth = 0.48),
axis.ticks.y = element_line(colour = "#334650", linewidth = 0.42),
axis.ticks.x = element_line(colour = "#334650", linewidth = 0.42),
axis.text.y = element_text(colour = "#263B47", face = "bold", size = 10.5,
hjust = 0.5, lineheight = 0.9,
margin = margin(r = 8)),
axis.text.x = element_text(colour = "#425763", size = 9.5),
axis.title.x = element_text(colour = "#263B47", face = "bold", size = 10.5,
margin = margin(t = 10)),
plot.caption = element_text(colour = "#5A6D77", size = 8.3,
hjust = 0, lineheight = 0.96,
margin = margin(t = 9)),
legend.position = "inside",
legend.position.inside = c(0.82, 0.065),
legend.justification = c(0.5, 0),
legend.background = element_rect(fill = "white", colour = NA),
legend.title = element_text(colour = "#263B47", face = "bold", size = 9.0),
legend.text = element_text(colour = "#425763", size = 8.3),
plot.margin = margin(15, 14, 12, 12)
)
raincloud_plot进阶图怎么读?
每个方形点是一个国家;半密度轮廓显示组内分布形状,其高度不表示该洲国家总数或人口数。红菱形是国家等权的均值,箱体中线是中位数。先对四洲做
Kruskal–Wallis 总体检验,再用 Dunn 检验比较各洲;括号显示经 Holm
校正后达到显著水平的五组比较。*、**、***
分别表示校正后 P <
0.05、0.01、0.001;亚洲与美洲未达显著水平。检验以国家为单位,不代表居民个人层面的差异。
1 · 图形与适用问题
流向图(alluvial plot)把同一批记录在多个分类变量中的归属并排显示:每一列是一个分类变量,矩形是类别,连接带代表同时属于这些类别的记录。带和矩形的纵向厚度对应人数或其他权重。若变量确实按时间排列,它也能展示状态转移;若变量只是不同属性,横向顺序仅是阅读安排。
拿到“多个分类变量 +
每种组合的人数或权重”的数据时,可用流向图观察组成、分叉与汇合。基础图按“舱位或船员身份—性别—生还结果”排列三个变量;进阶图再加入年龄组。横向顺序只用于阅读,不表示事件先后。本节使用
ggalluvial 的 geom_alluvium()
画连接带、geom_stratum() 画类别矩形。
2 · 数据处理
R 内置的 Titanic
是四维汇总频数表,Freq
是每组人数;Class 除了一、二、三等舱,还包含
Crew(船员),所以数据对象不能统称为“乘客”。先转换为数据框,基础图即可使用。
事件背景:
泰坦尼克号(RMS Titanic)是一艘远洋客轮。1912 年 4 月 14
日深夜在北大西洋撞上冰山,15 日凌晨沉没,造成大量人员伤亡。
数据是什么:
R 内置的 Titanic
是按乘客与船员类别组合统计的汇总频数表(contingency table),不是逐人记录。
数据包含什么:
数据包含
Class(舱位或船员身份)、Sex(性别)、Age(成人/儿童)、Survived(是否幸存);Freq
表示该组合的人数。
3 · 最小代码与基础图
数据已经准备好后,axis1、axis2、axis3
指定并列的分类列,y = Freq
指定人数。两种几何对象分别画连接带和节点;基础图暂不加入配色、标签与版式设置。
ggplot(titanic_df, aes(axis1 = Class, axis2 = Sex,
axis3 = Survived, y = Freq)) +
geom_alluvium() +
geom_stratum()泰坦尼克号乘客及船员按身份、性别和生还结果排列的基础流向图
先看节点,再沿带追踪组合。 左侧 Class
节点的纵向厚度是各类总人数;中间和右侧同理。某条带在纵向有多厚,就表示该类别组合有多少人;横向长短与人数无关。右侧
No 节点厚于
Yes,表示这份乘客与船员合计数据中遇难者更多。
本图展示的是人数,不同类别的基数不一样。若要比较各舱位或不同性别的生还率,需另算“生还人数
÷ 该组总人数”,不能只比带的厚度。Class → Sex → Survived
只是三个属性的并列排列,不是时间上的状态转换。
4 · 进阶成图与美化思路
进阶图在基础图中加入年龄组,按舱位/身份—性别—年龄组—结局排列四列。人数大于 0 的组合共有 24 条路径、2,201 人。图上的三个视觉编码分别承担不同任务:
读图时从左侧身份沿连接带追踪到结局。纵轴提供各列的累计人数刻度,比较单条路径的人数应看带的厚度。路径占比不是某个舱位、性别或年龄组内部的生还率;Age
只区分成人与儿童,不表示具体年龄。
泰坦尼克号乘客与船员按舱位或身份、性别、年龄组和结局排列的四轴流向图;连接带的纵向厚度表示人数,颜色区分舱位与船员身份,透明度表示完整路径占全体 2,201 人的比例。
# 进阶流向图:泰坦尼克号的四维分类路径
library(dplyr)
library(ggplot2)
library(ggalluvial)
# 1. 选用本机中文字体,不联网下载
font_priority <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Hiragino Sans GB", "Heiti SC"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
available <- unique(systemfonts::system_fonts()$family)
candidates <- font_priority[font_priority %in% available]
has_chinese <- function(family) {
tryCatch(
all(systemfonts::glyph_info(c("中", "文", "R"), family = family)$index > 0),
error = function(e) FALSE
)
}
candidates <- candidates[vapply(candidates, has_chinese, logical(1))]
if (length(candidates)) font_family <- candidates[[1]]
}
if (.Platform$OS.type == "windows" && font_family != "sans") {
try(do.call(
grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(font_family)), font_family)
), silent = TRUE)
}
# 2. 每行是一条「身份 × 性别 × 年龄组 × 结局」路径
titanic_adv <- as.data.frame(Titanic) %>%
filter(Freq > 0) %>%
mutate(
prop = Freq / sum(Freq),
Class = factor(Class, levels = c("1st", "2nd", "3rd", "Crew"),
labels = c("一等", "二等", "三等", "船员")),
Sex = factor(Sex, levels = c("Female", "Male"),
labels = c("女性", "男性")),
Age = factor(Age, levels = c("Child", "Adult"),
labels = c("儿童", "成人")),
Survived = factor(Survived, levels = c("No", "Yes"),
labels = c("遇难", "生还"))
)
stopifnot(nrow(titanic_adv) == 24L, sum(titanic_adv$Freq) == 2201L)
# 3. 保留旧图的 NPG 四色,透明度另表示单条路径所占比例
class_cols <- ggsci::pal_npg()(4)
names(class_cols) <- levels(titanic_adv$Class)
plot_sankey <- ggplot(
titanic_adv,
aes(axis1 = Class, axis2 = Sex, axis3 = Age,
axis4 = Survived, y = Freq)
) +
geom_alluvium(
aes(fill = Class, alpha = prop),
width = 1 / 8, colour = NA
) +
geom_stratum(
width = 1 / 8,
fill = "#F5F7F8", colour = "#66747E", linewidth = 0.38
) +
geom_text(
stat = "stratum",
aes(label = after_stat(stratum)),
family = font_family, size = 3.35, colour = "#26343C"
) +
scale_x_discrete(
limits = c("Class", "Sex", "Age", "Survived"),
labels = c(Class = "舱位/身份", Sex = "性别",
Age = "年龄组", Survived = "结局"),
expand = expansion(add = c(0.18, 0.18))
) +
scale_y_continuous(
breaks = seq(0, 2000, by = 500),
labels = scales::label_comma(),
expand = expansion(mult = c(0, 0.015))
) +
scale_fill_manual(name = "舱位/身份", values = class_cols) +
scale_alpha_continuous(
name = "路径占比(透明度)",
limits = c(0, 0.31),
breaks = c(0.01, 0.05, 0.10, 0.20, 0.30),
labels = scales::label_percent(accuracy = 1),
range = c(0.22, 0.78)
) +
guides(
fill = guide_legend(
order = 1,
override.aes = list(alpha = 0.78, colour = NA)
),
alpha = guide_legend(
order = 2,
override.aes = list(fill = "#3D667A", colour = NA)
)
) +
labs(x = NULL, y = "累计人数(人)") +
theme_classic(base_family = font_family, base_size = 10) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
panel.grid.major.y = element_line(colour = "#E9EEF1", linewidth = 0.22),
axis.line = element_line(colour = "#596873", linewidth = 0.38),
axis.ticks = element_line(colour = "#596873", linewidth = 0.33),
axis.ticks.length = grid::unit(1.8, "mm"),
axis.title.y = element_text(colour = "#2B3B46", size = 9.1,
margin = margin(r = 8)),
axis.text.y = element_text(colour = "#4F606B", size = 8.4),
axis.text.x = element_text(colour = "#2B3B46", face = "bold",
size = 9.2, margin = margin(t = 5)),
legend.position = "right",
legend.justification = "top",
legend.box = "vertical",
legend.box.spacing = grid::unit(4, "mm"),
legend.spacing.y = grid::unit(1.5, "mm"),
legend.title = element_text(colour = "#2B3B46", face = "bold", size = 9),
legend.text = element_text(colour = "#40515D", size = 8.7),
legend.key.width = grid::unit(11, "mm"),
legend.key.height = grid::unit(4.4, "mm"),
legend.margin = margin(0, 0, 0, 0),
plot.margin = margin(8, 10, 8, 5)
)
print(plot_sankey)1 · 图形与适用问题
山峦图(ridgeline plot)把多个组别的连续变量分布画成上下错开的密度曲线。手中有“一个连续变量 + 一个分组变量”,且希望比较多组的典型取值、分布宽窄或峰形时,可以考虑它。若分组是不同年份,也可用来观察分布随时间怎样移动。
每一条“山脊”是核密度估计(KDE):算法将样本点平滑成一条曲线。曲线形状受样本量和带宽影响,不能把某个小突起直接当成稳定的人群分层;山体之间的视觉重叠也不能单独证明统计差异“显著”。
2 · 数据处理
本例选取 gapminder 的 2007
年数据,以国家为观测单位,比较各洲国家的预期寿命。大洋洲在该年只有两个国家,不足以稳定估计密度,因此本例不画其密度曲线;这是样本限制,不是大洋洲没有数据。
library(ggplot2)
library(gapminder)
library(ggridges)
ridge_data <- subset(gapminder, year == 2007 & continent != "Oceania")3 · 最小代码与基础图
只需映射连续的 lifeExp 与分组的
continent,再加入 geom_density_ridges()。
2007 年四大洲国家预期寿命分布的基础山峦图
本图中,欧洲的密度主要位于较高的寿命区间,非洲整体偏向较低区间;亚洲和美洲有较宽的取值范围。各洲国家数不同,曲线展示的是组内分布形状,山体高度不能直接当作该洲人口或国家总数。
4 · 进阶成图与美化思路
进阶图聚焦非洲与亚洲,使用 1952—2007 年每五年一次的国家记录。每条山脊对应某一大洲、某一年的国家等权预期寿命分布;两列共用横轴、密度高度尺度和 3.4 岁的核密度带宽。大洋洲每年只有 2 个国家,不宜估计平滑密度;美洲和欧洲未纳入这次双洲对照。
读图时,先沿年份比较分布的位置与宽窄,再用山脊下方的短刻度核对各国观测。珊瑚色点和旁边的数字表示当年的国家中位数,虚线连接这些中位数。填色采用 viridisLite 的 Mako 色盘:预期寿命越高颜色越深,并在较高寿命区间分配更密的色阶;颜色差不能直接当作岁数差。相邻山脊的透明重叠不表示国家数或显著性,曲线高度也不能当作人口或国家总数。
1952—2007 年非洲与亚洲各国出生时预期寿命的分布。每五年一条山脊;短刻度为各国观测,珊瑚色点和数字为逐年中位数。密度带宽统一为 3.4 岁;大洋洲每年仅 2 国,未纳入本图。
# 进阶山峦图:非洲与亚洲各国预期寿命的逐年分布
suppressPackageStartupMessages({
library(dplyr)
library(ggplot2)
library(gapminder)
library(tidyr)
library(purrr)
})
# 1. 字体 ---------------------------------------------------------------
font_priority <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Hiragino Sans GB", "Heiti SC"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
installed <- unique(systemfonts::system_fonts()$family)
candidates <- font_priority[font_priority %in% installed]
has_chinese <- function(family) {
tryCatch(
all(
systemfonts::glyph_info(
c("中", "文", "R"),
family = family
)$index > 0
),
error = function(e) FALSE
)
}
candidates <- candidates[
vapply(candidates, has_chinese, logical(1))
]
if (length(candidates)) {
font_family <- candidates[[1]]
}
}
if (.Platform$OS.type == "windows" && font_family != "sans") {
try(
do.call(
grDevices::windowsFonts,
setNames(
list(grDevices::windowsFont(font_family)),
font_family
)
),
silent = TRUE
)
}
# 2. 数据 ---------------------------------------------------------------
years <- sort(unique(gapminder$year))
ridge_data <- gapminder %>%
filter(continent %in% c("Africa", "Asia")) %>%
mutate(
continent = factor(
as.character(continent),
levels = c("Africa", "Asia")
),
year_id = match(year, years)
) %>%
arrange(continent, year, country)
ridge_stats <- ridge_data %>%
group_by(continent, year, year_id) %>%
summarise(
median_life = median(lifeExp),
n_country = n_distinct(country),
.groups = "drop"
) %>%
arrange(continent, year_id)
stopifnot(
length(years) == 12L,
nrow(ridge_data) == 1020L,
all(ridge_stats$n_country[ridge_stats$continent == "Africa"] == 52L),
all(ridge_stats$n_country[ridge_stats$continent == "Asia"] == 33L),
!anyNA(ridge_data$lifeExp)
)
# 3. 参数 ---------------------------------------------------------------
x_min <- 18
x_max <- 90
density_bw <- 3.4
n_grid <- 900
ridge_overlap <- 1.60
fill_alpha <- 0.82
# 色块边界略微交叠,避免显示时出现细缝
x_bleed_mult <- 0.55 # 横向小片轻微重叠
top_bleed <- 0.025 # 略微覆盖轮廓边缘,避免渲染细缝
# 反向 Mako 色盘:预期寿命越高,颜色越深
# 色标节点并非等距,较高寿命区间有更多颜色变化
fill_knots <- c(18, 40, 50, 58, 65, 70, 75, 80, 90)
fill_cols <- scales::alpha(
viridisLite::mako(9, begin = 0.08, end = 0.92, direction = -1),
fill_alpha
)
# 4. KDE 计算 -----------------------------------------------------------
ridge_stats <- ridge_stats %>%
mutate(
median_label = sprintf("%.1f", median_life),
median_label_colour = if_else(median_life >= 64, "#FFFFFF", "#203846")
)
ridge_density <- ridge_data %>%
group_by(continent, year, year_id) %>%
group_modify(~ {
d <- density(
.x$lifeExp,
bw = density_bw,
from = x_min,
to = x_max,
n = n_grid
)
tibble(
x = d$x,
dens = d$y
)
}) %>%
ungroup()
# 共同高度尺度
height_mult <- ridge_overlap / max(ridge_density$dens)
ridge_density <- ridge_density %>%
mutate(
y0 = year_id,
y = year_id + dens * height_mult
)
# 5. 构造连续渐变填充面 -------------------------------------------------
# 小色块在横向和顶部略微相接,形成连续填色
ridge_fill <- ridge_density %>%
group_by(continent, year, year_id) %>%
group_modify(~ {
d <- arrange(.x, x)
x1 <- d$x[-nrow(d)]
x2 <- d$x[-1]
y1 <- d$y[-nrow(d)]
y2 <- d$y[-1]
y0 <- d$y0[1]
dx <- x2 - x1
x1_bleed <- pmax(x_min, x1 - dx * x_bleed_mult)
x2_bleed <- pmin(x_max, x2 + dx * x_bleed_mult)
# 顶部填色略微相接
y1_bleed <- y1 + top_bleed
y2_bleed <- y2 + top_bleed
tibble(
piece_id = paste0(as.character(.y$continent), "_", .y$year, "_", seq_along(x1)),
fill_x = (x1 + x2) / 2,
x = pmap(
list(x1_bleed, x2_bleed, x2_bleed, x1_bleed),
c
),
y = pmap(
list(y0, y0, y2_bleed, y1_bleed),
c
)
) %>%
unnest(c(x, y))
}) %>%
ungroup()
# 6. 作图 ---------------------------------------------------------------
plot_ridge <- ggplot() +
# 背景横线
geom_hline(
yintercept = seq_along(years),
colour = "#E7EDF1",
linewidth = 0.22
) +
# 连续渐变填充面(先画 fill)
geom_polygon(
data = ridge_fill,
aes(
x = x,
y = y,
group = piece_id,
fill = fill_x
),
colour = NA
) +
# 真实观测值
geom_segment(
data = ridge_data,
aes(
x = lifeExp,
xend = lifeExp,
y = year_id - 0.10,
yend = year_id - 0.025
),
colour = "#3E5562",
linewidth = 0.26,
alpha = 0.38
) +
# 中位数轨迹
geom_path(
data = ridge_stats,
aes(
x = median_life,
y = year_id,
group = continent,
linetype = "中位数"
),
colour = "#253B47",
linewidth = 0.45,
lineend = "round"
) +
# 中位数点
geom_point(
data = ridge_stats,
aes(
x = median_life,
y = year_id,
shape = "中位数"
),
fill = "#F16C4A",
colour = "white",
size = 2.2,
stroke = 0.55
) +
# 山脊轮廓置于填色上方
geom_line(
data = ridge_density,
aes(
x = x,
y = y,
group = interaction(continent, year)
),
colour = "#314956",
linewidth = 0.34,
lineend = "round"
) +
# 每个年度中位数的具体岁数;深色区域使用白字
geom_text(
data = ridge_stats,
aes(
x = median_life + 1.3,
y = year_id + 0.11,
label = median_label,
colour = median_label_colour
),
hjust = 0,
vjust = 0.5,
family = font_family,
fontface = "bold",
size = 2.5,
show.legend = FALSE
) +
facet_wrap(
~ continent,
nrow = 1,
labeller = as_labeller(
c(
Africa = "非洲 · 每年 52 国",
Asia = "亚洲 · 每年 33 国"
)
)
) +
scale_fill_gradientn(
colours = fill_cols,
values = scales::rescale(fill_knots, from = c(x_min, x_max)),
limits = c(x_min, x_max),
breaks = c(30, 50, 60, 70, 80),
name = "预期寿命(岁)"
) +
scale_colour_identity() +
scale_shape_manual(
values = c("中位数" = 21),
name = NULL
) +
scale_linetype_manual(
values = c("中位数" = "22"),
name = NULL
) +
scale_x_continuous(
breaks = seq(20, 90, by = 10),
limits = c(x_min, x_max),
expand = expansion(mult = c(0, 0))
) +
scale_y_continuous(
breaks = seq_along(years),
labels = years,
limits = c(0.65, 14.25),
expand = expansion(mult = c(0, 0))
) +
guides(
fill = guide_colourbar(order = 1, barheight = grid::unit(42, "mm")),
shape = guide_legend(order = 2),
linetype = guide_legend(order = 2)
) +
labs(
x = "出生时预期寿命(岁)",
y = NULL
) +
theme_classic(
base_size = 10,
base_family = font_family
) +
theme(
plot.background = element_rect(
fill = "white",
colour = NA
),
panel.background = element_rect(
fill = "white",
colour = NA
),
panel.grid.major.x = element_line(
colour = "#E9EEF1",
linewidth = 0.22
),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.line.x = element_line(
colour = "#73848E",
linewidth = 0.30
),
axis.line.y = element_blank(),
axis.ticks.x = element_line(
colour = "#73848E",
linewidth = 0.27
),
axis.ticks.y = element_blank(),
axis.text = element_text(
colour = "#344C5B",
size = 8.8
),
axis.title.x = element_text(
colour = "#253D4B",
face = "bold",
size = 10,
margin = margin(t = 7)
),
strip.background = element_blank(),
strip.text = element_text(
colour = "#253D4B",
face = "bold",
size = 10,
margin = margin(t = 2, b = 7)
),
panel.spacing = grid::unit(0.75, "lines"),
legend.position = "right",
legend.title = element_text(
colour = "#253D4B",
face = "bold",
size = 9.2
),
legend.text = element_text(
colour = "#425663",
size = 8.7
),
legend.key.width = grid::unit(0.9, "cm"),
plot.margin = margin(8, 10, 7, 8)
)
stopifnot(dplyr::n_distinct(ridge_fill$piece_id) == nrow(ridge_stats) * (n_grid - 1L))
print(plot_ridge)1 · 图形与适用问题
热力图(heatmap)把两个维度的组合排成网格,以颜色表示每格的数值。手中有“行分类 × 列分类(或时间) × 数值”的数据时,适合用它概览高低区域、长期变化和分组差异;把多个数值变量两两配对,还可以用热力图展示相关系数。颜色利于快速发现模式;要读出精确数值,仍需查看格内标注、色标或原始数据。
本例横轴是年份,纵轴是大洲,颜色代表该年该洲的人口加权平均预期寿命。同一套色标贯穿整张图,才能比较不同行与列。
2 · 数据处理
gapminder 的每一行是一个国家在某一年的记录。绘制“大洲 ×
年份”热力图前,先按这两个变量汇总。人口加权均值的计算方法为
sum(lifeExp * pop) / sum(pop);它让人口较多的国家在洲级指标中占更大权重。这里的世界和大洲统计仅覆盖
gapminder 收录的国家与年份。
library(dplyr)
library(ggplot2)
library(gapminder)
heat_df <- gapminder %>%
group_by(continent, year) %>%
summarise(life_exp_w = weighted.mean(lifeExp, pop),
.groups = "drop")3 · 最小代码与基础图
数据已有“每洲每年一行”,只需映射网格位置和填充值,再调用
geom_tile()。基础图先保留默认色标。
各洲历年人口加权平均预期寿命的基础热力图(单位:岁)
若想判断某段时间是否“加速”或“放缓”,要回到数值计算相邻年份的变化率;单凭渐变颜色不够。图中每格是人口加权均值,不代表该洲所有国家都处在同一水平。
4 · 进阶成图与美化思路
下面用两张图展示热力图的不同用途:第一张沿日期查找高低值,第二张比较变量之间的关联。
进阶示例 1 / 2 · 日历热力图
基础图按“大洲 × 年份”比较汇总值。进阶图改用真正的日历布局:2023 年的 12 个月各占一个面板,横向为周一至周日,纵向为月内周次,每格代表一天。颜色表示当天有起飞延误记录的航班中,起飞延误至少 15 分钟的比例。所有月份共用同一套百分比色标。
两张进阶图使用公开 R 包
nycflights23:flights 每行是一趟 2023
年纽约地区机场的航班记录,weather
提供机场天气记录。首次使用时,请在 Console 运行
install.packages(c("nycflights23", "viridisLite", "corrplot", "colorspace"));安装是一次性操作,不写入绘图代码。查看 CRAN
数据包说明。
dep_delay
缺失的航班不算作准点航班,也不进入本图每日比例的分母。
先准备数据。
按日期计算延误航班数与有起飞延误记录的航班数,再用
seq.Date() 列出全年 365
天。计算“星期列”和“月内周次”后,数据中每一天才会落在正确的日历位置。
再看成图。 色块越深,当天达到延误阈值的航班比例越高;颜色表示的是比例,不代表延误航班的总数,也不能单凭图形判断延误原因。
2023 年纽约地区机场出发航班逐日起飞延误至少 15 分钟的比例。每格代表一天;分母为当日有起飞延误记录的航班。
先准备数据:计算每日比例与日历位置。
if (.Platform$OS.type == "windows") {
invisible(try(Sys.setlocale("LC_CTYPE", ".UTF-8"), silent = TRUE))
}
library(dplyr)
library(ggplot2)
# 1. 每天一行:分母只包含有起飞延误记录的航班。
calendar_daily <- nycflights23::flights %>%
filter(year == 2023, !is.na(dep_delay)) %>%
mutate(
date = as.Date(sprintf("%04d-%02d-%02d", year, month, day))
) %>%
group_by(date) %>%
summarise(
n_departed = n(),
n_delayed = sum(dep_delay >= 15),
delay_rate = n_delayed / n_departed,
.groups = "drop"
)
# 2. 先列出全年日期,再合并每日比例;无记录日期不会悄悄消失。
calendar_df <- tibble(
date = seq.Date(
as.Date("2023-01-01"),
as.Date("2023-12-31"),
by = "day"
)
) %>%
left_join(calendar_daily, by = "date") %>%
mutate(
month = factor(
as.integer(format(date, "%m")),
levels = 1:12,
labels = paste0(1:12, "月")
),
day_number = as.integer(format(date, "%d")),
# R 默认周日=0;转换成周一=1、周日=7。
weekday = (as.POSIXlt(date)$wday + 6L) %% 7L + 1L,
month_start = as.Date(format(date, "%Y-%m-01")),
first_weekday =
(as.POSIXlt(month_start)$wday + 6L) %% 7L + 1L,
# 月初按星期向右错位;每经过七天转到下一行。
week_in_month =
(day_number + first_weekday - 2L) %/% 7L + 1L,
calendar_y = 7L - week_in_month
)
stopifnot(
nrow(calendar_df) == 365L,
!anyNA(calendar_df$delay_rate),
all(calendar_df$n_departed > 0)
)
weekday_header <- tibble(
weekday = 1:7,
calendar_y = 7,
label = c(
"周一", "周二", "周三", "周四",
"周五", "周六", "周日"
)
)再绘图:设置日期格、颜色与月份布局。
# 1. 字体:优先使用本机中文无衬线字体,不联网下载。
font_priority <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Hiragino Sans GB", "Heiti SC"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
installed <- unique(systemfonts::system_fonts()$family)
candidates <- font_priority[font_priority %in% installed]
has_chinese <- function(family) {
tryCatch(
all(systemfonts::glyph_info(
c("月", "日", "延"), family = family
)$index > 0),
error = function(e) FALSE
)
}
candidates <- candidates[vapply(candidates, has_chinese, logical(1))]
if (length(candidates)) font_family <- candidates[[1]]
}
if (.Platform$OS.type == "windows" && font_family != "sans") {
try(
do.call(
grDevices::windowsFonts,
setNames(
list(grDevices::windowsFont(font_family)),
font_family
)
),
silent = TRUE
)
}
# 2. 顺序色标与日期字色:比例越高,色块越深。
colour_max <- ceiling(max(calendar_df$delay_rate) * 10) / 10
fill_colours <- viridisLite::mako(
9, begin = 0.08, end = 0.96, direction = -1
)
calendar_df <- calendar_df %>%
mutate(
label_colour = if_else(
delay_rate >= 0.36, "#FFFFFF", "#000000"
)
)
# 3. 每个月复用同一套星期位置和百分比色标。
plot_calendar <- ggplot(calendar_df, aes(weekday, calendar_y)) +
geom_tile(
aes(fill = delay_rate),
width = 0.94, height = 0.94,
colour = "#8499A2", linewidth = 0.26
) +
geom_text(
aes(label = day_number, colour = label_colour),
family = font_family, size = 2.55,
show.legend = FALSE
) +
geom_text(
data = weekday_header,
aes(x = weekday, y = calendar_y, label = label),
inherit.aes = FALSE,
family = font_family, fontface = "bold",
size = 2.25, colour = "#5C6D76"
) +
facet_wrap(~ month, ncol = 4, drop = FALSE) +
scale_fill_gradientn(
colours = fill_colours,
limits = c(0, colour_max),
breaks = seq(0, colour_max, by = 0.1),
labels = scales::label_percent(accuracy = 1),
name = "延误 ≥15 分钟"
) +
scale_colour_identity() +
scale_x_continuous(
limits = c(0.5, 7.5), expand = c(0, 0)
) +
scale_y_continuous(
limits = c(0.5, 7.5), expand = c(0, 0)
) +
coord_fixed(clip = "off") +
guides(
fill = guide_colourbar(
title.position = "top",
barheight = grid::unit(46, "mm"),
barwidth = grid::unit(4.3, "mm")
)
) +
labs(title = "2023", x = NULL, y = NULL) +
theme_void(base_family = font_family, base_size = 10) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA),
plot.title.position = "plot",
plot.title = element_text(
hjust = 0.5, face = "bold", size = 15,
colour = "#253F4C", margin = margin(b = 12)
),
strip.background = element_blank(),
strip.text = element_text(
colour = "#253F4C", face = "bold", size = 10.2,
margin = margin(b = 5)
),
panel.spacing.x = grid::unit(0.9, "lines"),
panel.spacing.y = grid::unit(0.6, "lines"),
legend.position = "right",
legend.box.spacing = grid::unit(9, "mm"),
legend.title = element_text(
colour = "#253F4C", face = "bold", size = 9
),
legend.text = element_text(
colour = "#465B67", size = 8.5
),
plot.margin = margin(8, 10, 7, 8)
)
print(plot_calendar)进阶示例 2 / 2 · 相关性热力图
热力图还可以展示多个数值变量的两两关系。上一张日历图覆盖纽约地区三个出发机场;这一张聚焦 JFK 机场,把 2023 年的每日起飞延误比例与六项天气指标按日期对齐。七个指标形成 21 组不重复的比较,因此图中只画下三角,也不显示变量与自身的相关。
ns
对应校正后的检验结果。星号越多不等于关联越强,更不能据此认定因果关系。先准备数据。 flights
计算每日延误比例,weather 汇总同日天气;天气数据未覆盖 12
月 31 日,所以最终使用 364
个机场日。代码先将变量转为秩,再扣除年内趋势与星期效应,避免把这些共同的时间变化直接当作变量之间的关系。
再看成图。
紫色表示负相关,绿色表示正相关;格内数字给出调整后的相关系数,数字下方的星号或
ns 表示检验结果。
JFK 机场 2023 年每日延误比例与天气指标的调整后秩相关系数。只显示不重复的下三角比较。
读图示例。 延误比例与降水量的系数约为
0.49,与能见度的系数约为
−0.45;前者为同向关联,后者为反向关联。气温与延误比例接近 0,标为
ns。这些是同日数据的关联,不能据此确定延误原因。
先准备数据:汇总机场日并计算相关系数与检验结果。
library(dplyr)
library(ggplot2)
if (.Platform$OS.type == "windows") {
invisible(try(Sys.setlocale("LC_CTYPE", ".UTF-8"), silent = TRUE))
}
# 1. 按日汇总;延误定义与前面的日历图一致:起飞延误至少 15 分钟。
mean_or_na <- function(x) if (all(is.na(x))) NA_real_ else mean(x, na.rm = TRUE)
sum_or_na <- function(x) if (all(is.na(x))) NA_real_ else sum(x, na.rm = TRUE)
flight_daily <- nycflights23::flights %>%
filter(origin == "JFK", year == 2023, !is.na(dep_delay)) %>%
mutate(date = as.Date(sprintf("%04d-%02d-%02d", year, month, day))) %>%
group_by(date) %>%
summarise(delay_rate = mean(dep_delay >= 15), .groups = "drop")
weather_daily <- nycflights23::weather %>%
filter(origin == "JFK", year == 2023) %>%
mutate(date = as.Date(sprintf("%04d-%02d-%02d", year, month, day))) %>%
group_by(date) %>%
summarise(
temp_c = (mean_or_na(temp) - 32) * 5 / 9,
humid_pct = mean_or_na(humid),
wind_kmh = mean_or_na(wind_speed) * 1.60934,
rain_mm = sum_or_na(precip) * 25.4,
pressure_hpa = mean_or_na(pressure),
visib_km = mean_or_na(visib) * 1.60934,
.groups = "drop"
)
corr_daily <- inner_join(flight_daily, weather_daily, by = "date") %>%
arrange(date) %>%
filter(if_all(-date, is.finite))
stopifnot(nrow(corr_daily) == 364L) # 天气数据没有 2023-12-31
corr_input <- corr_daily %>%
select(delay_rate, rain_mm, humid_pct, visib_km,
pressure_hpa, temp_c, wind_kmh)
names(corr_input) <- c(
"延误比例", "降水量", "相对湿度", "能见度", "气压", "气温", "风速"
)
# 2. 先取秩,再扣除年内平滑趋势与星期效应。
# 得到的矩阵反映“同一天的相对变化”之间的关联。
rank_mat <- apply(corr_input, 2, rank, ties.method = "average")
day_index <- as.integer(corr_daily$date - min(corr_daily$date)) + 1L
weekday <- factor(as.integer(format(corr_daily$date, "%u")), levels = 1:7)
basis <- model.matrix(~ splines::ns(day_index, df = 6) + weekday)
z <- scale(lm.fit(basis, rank_mat)$residuals)
stopifnot(all(is.finite(z)))
corr_mat <- crossprod(z) / (nrow(z) - 1)
dimnames(corr_mat) <- list(names(corr_input), names(corr_input))
# 3. 显著性:连续 28 天为一区块,只重排区块顺序,保留区块内日期顺序。
# 双侧置换检验;对所有 21 组配对统一做 Holm 多重比较校正。
pair_id <- which(lower.tri(corr_mat), arr.ind = TRUE)
rho_obs <- corr_mat[pair_id]
blocks <- split(seq_len(nrow(z)), ceiling(seq_len(nrow(z)) / 28L))
n_perm <- 99999L
extreme <- integer(nrow(pair_id))
set.seed(927)
for (b in seq_len(n_perm)) {
new_order <- unlist(blocks[sample.int(length(blocks))], use.names = FALSE)
perm_corr <- crossprod(z, z[new_order, , drop = FALSE]) / (nrow(z) - 1)
extreme <- extreme + (abs(perm_corr[pair_id]) >= abs(rho_obs))
}
p_raw <- (extreme + 1) / (n_perm + 1)
p_adj <- p.adjust(p_raw, method = "holm")
p_mat <- matrix(NA_real_, nrow(corr_mat), ncol(corr_mat),
dimnames = dimnames(corr_mat))
p_mat[pair_id] <- p_adj
p_mat[cbind(pair_id[, 2], pair_id[, 1])] <- p_adj
# 聚类只决定变量显示顺序,不改动相关系数或显著性结果。
order_id <- corrplot::corrMatOrder(
corr_mat, order = "hclust", hclust.method = "average"
)
var_order <- colnames(corr_mat)[order_id]
corr_long <- as.data.frame(as.table(corr_mat), stringsAsFactors = FALSE)
names(corr_long) <- c("row_var", "col_var", "rho")
corr_long <- corr_long %>%
mutate(
row_var = as.character(row_var),
col_var = as.character(col_var),
row_id = match(row_var, var_order),
col_id = match(col_var, var_order),
p_adj = p_mat[cbind(row_var, col_var)]
) %>%
filter(row_id > col_id) %>%
mutate(
stars = case_when(
p_adj < 0.001 ~ "***",
p_adj < 0.01 ~ "**",
p_adj < 0.05 ~ "*",
TRUE ~ "ns"
),
number = sprintf("%.2f", rho),
col_var = factor(col_var, levels = var_order[-length(var_order)]),
row_var = factor(row_var, levels = rev(var_order[-1]))
)
stopifnot(nrow(corr_long) == 21L, !anyNA(corr_long$p_adj))再绘图:设置下三角、发散色标与格内标注。
# 4. 绿色—紫色发散色盘:负值紫、零值浅灰、正值绿。
pal_corr <- colorspace::diverging_hcl(256, palette = "Purple-Green")
tile_colour <- scales::gradient_n_pal(pal_corr)((corr_long$rho + 1) / 2)
contrast_white <- colorspace::contrast_ratio("#FFFFFF", tile_colour)
contrast_black <- colorspace::contrast_ratio("#000000", tile_colour)
corr_long$text_colour <- ifelse(
contrast_white >= contrast_black, "#FFFFFF", "#000000"
)
stopifnot(min(pmax(contrast_white, contrast_black)) >= 4.5)
# 优先使用系统自带的中文无衬线字体;不联网下载字体。
font_priority <- switch(
Sys.info()[["sysname"]],
Windows = c("Microsoft YaHei", "Microsoft YaHei UI", "SimHei"),
Darwin = c("PingFang SC", "Hiragino Sans GB", "Heiti SC"),
Linux = c("Noto Sans CJK SC", "WenQuanYi Micro Hei"),
"sans"
)
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
installed <- unique(systemfonts::system_fonts()$family)
for (candidate in font_priority[font_priority %in% installed]) {
has_chinese <- tryCatch(
all(systemfonts::glyph_info(c("延", "误", "降"), family = candidate)$index > 0),
error = function(e) FALSE
)
if (has_chinese) {
font_family <- candidate
break
}
}
}
if (.Platform$OS.type == "windows" && font_family != "sans") {
try(do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(font_family)), font_family)),
silent = TRUE)
}
# 5. 成图:下三角、底部列名、数字和检验标注分行。
plot_corr <- ggplot(corr_long, aes(x = col_var, y = row_var)) +
geom_tile(
aes(fill = rho),
width = 0.84, height = 0.84,
colour = NA
) +
geom_text(
aes(label = number, colour = text_colour),
position = position_nudge(y = 0.14),
family = font_family, fontface = "bold", size = 3.4,
show.legend = FALSE
) +
geom_text(
data = filter(corr_long, stars != "ns"),
aes(label = stars, colour = text_colour),
position = position_nudge(y = -0.18),
family = font_family, fontface = "bold", size = 3.0,
show.legend = FALSE
) +
geom_text(
data = filter(corr_long, stars == "ns"),
aes(label = stars, colour = text_colour),
position = position_nudge(y = -0.18),
family = font_family, fontface = "plain", size = 2.85,
show.legend = FALSE
) +
scale_colour_identity(guide = "none") +
scale_fill_gradientn(
colours = pal_corr,
limits = c(-1, 1),
breaks = c(-1, -0.5, 0, 0.5, 1),
name = "调整后秩相关系数",
guide = guide_colourbar(
title.position = "top",
barheight = grid::unit(43, "mm"),
barwidth = grid::unit(4.5, "mm"),
frame.colour = NA,
ticks.colour = "#506168"
)
) +
scale_x_discrete(position = "bottom", expand = expansion(add = 0.24)) +
scale_y_discrete(expand = expansion(add = 0.24)) +
coord_fixed(clip = "off") +
labs(
x = NULL, y = NULL,
caption = paste0(
"注:n = 364 个机场日。各变量先取秩,扣除年内趋势和星期效应后计算残差相关;同期相关不代表因果。\n",
"* 校正后 p < 0.05;** < 0.01;*** < 0.001;ns:校正后 p ≥ 0.05。"
)
) +
theme_minimal(base_family = font_family, base_size = 10.5) +
theme(
plot.background = element_rect(fill = "#FFFFFF", colour = NA),
panel.background = element_rect(fill = "#FFFFFF", colour = NA),
panel.grid = element_blank(),
axis.text.x = element_text(
colour = "#293E45", face = "bold", size = 9.6,
margin = margin(t = 8)
),
axis.text.y = element_text(
colour = "#293E45", face = "bold", size = 9.6,
margin = margin(r = 8)
),
axis.ticks = element_blank(),
legend.position = "right",
legend.title = element_text(
colour = "#263B41", face = "bold", size = 9.2
),
legend.text = element_text(colour = "#455860", size = 8.6),
plot.caption.position = "plot",
plot.caption = element_text(
colour = "#53636A", size = 8.0, hjust = 0,
lineheight = 1.22, margin = margin(t = 12)
),
plot.margin = margin(8, 10, 8, 8)
)
print(plot_corr)当分析重点集中在少数几个指标时,可以把它们与其他变量的相关系数集中展示。本例以售价和建成年份为核心指标,分别与同样的
16 项住房特征比较,便于查看同一特征在两个指标上的关联差异。图形通过
ggplot2 自定义圆周布局绘制。
数据来自美国爱荷华州 Ames 市 2006 至 2010 年的住宅交易资料,原表包含 2,930 套住宅。本例选取售价、建成年份及 16 项住房特征;排除所选字段存在缺失的 25 行后,统一使用 2,905 行完整记录。零值与合法极值保留,售价采用原表美元价格,未作通胀调整。
本例计算 Spearman 秩相关系数 ρ,描述两个变量数值排序上的单调关联;32 项检验统一进行 BH 多重比较校正,以减少同时检验多组关系时的误报。
Ames 住宅资料的环形相关性图:上半圈为售价,下半圈为建成年份,分别与同样的 16 项住房特征比较;圆内标出 Spearman ρ 与 BH 校正后的显著性星号。
按这个顺序读图。 点击图片可放大查看,窄屏可在图框内左右滑动。
| 阅读顺序 | 对应含义 |
|---|---|
| 先看半圈 | 上半圈对应售价,下半圈对应建成年份。 |
| 再找变量 | 同一住房特征上下镜像排列,外环短线帮助标签与圆点对应。 |
| 看数字与颜色 | 圆内数字为 ρ;红色为正相关,蓝色为负相关,接近白色表示接近 0。绝对值越大,单调关联越强。 |
| 最后看星号 | *、**、***
分别对应 BH 校正后 p < 0.05、0.01、0.001;没有星号表示未达到 0.05
阈值。 |
所有圆点大小相同,圆周位置和外环短线均不表示相关强弱。星号反映检验结果,关联强弱仍应看 ρ。 BH 校正只调整 p 值,不改变相关系数。
读图示例。 二层面积与售价的 ρ 约为 0.25,与建成年份约为 0.02,后者没有星号。比较上下对应的圆点,可以看出同一特征与两个核心指标的关联并不相同。
palette_colors
控制发散色盘,figure_width、figure_height
控制输出尺寸;更换颜色时保留以 0 为中心、范围为 −1 至 1 的色标。脚本将 PNG、PDF、相关系数表和变量说明写入 output
文件夹。整段复制运行时,data 和 output
建在当前工作目录;通过 source() 或 Rscript
运行时,建在脚本所在目录。年份和计数变量有重复取值,因此设置
exact = FALSE,使用近似检验。
以下安装命令在 Console 中运行一次;后面的脚本包含每次需要执行的包调用、数据处理、计算和导出步骤。
完整代码读取同一份作图脚本,末尾在交互会话中调用
print(),可在 RStudio
中查看图形。字体优先从本机已安装的中文字体中选择;如报“未找到可显示中文的字体”,先核对脚本的
font_candidates。
# Ames 房屋数据:单圆环形相关性图
# 运行方式:Rscript build-ames-circular-correlation.R
# 或在 RStudio 中打开本脚本,逐段运行;整段复制运行时,结果保存在当前工作目录。
# 首次安装仅在 Console 运行一次(本脚本不会自动安装包):
# install.packages(c("ggplot2", "dplyr", "cowplot", "ragg", "systemfonts", "jsonlite", "digest"))
# 环形排布借鉴公众号示例;数据与统计结果均重新计算。
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(cowplot)
})
# 1. 输出目录与可调整参数 -----------------------------------------------
get_script_dir <- function() {
source_files <- lapply(sys.frames(), function(frame) frame$ofile)
source_files <- Filter(function(x) !is.null(x) && nzchar(x), source_files)
if (length(source_files)) {
return(dirname(normalizePath(tail(source_files, 1)[[1]], mustWork = TRUE)))
}
file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
if (length(file_arg)) {
return(dirname(normalizePath(sub("^--file=", "", file_arg[1]), mustWork = TRUE)))
}
getwd()
}
script_dir <- get_script_dir()
data_dir <- file.path(script_dir, "data")
output_dir <- file.path(script_dir, "output")
dir.create(data_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
show_significance <- TRUE # 显示 BH 校正后的显著性星号
figure_width <- 10 # 英寸:一个圆,上下半圈分别对应两个核心指标
figure_height <- 10
figure_dpi <- 240
marker_center_radius <- 1.04 # 标记圆心所在的圆周
marker_radius <- 0.093 # 所有圆形标记等半径;大小不编码相关强弱
outer_ring_radius <- 1.19 # 外边界与标记分开,完整围住所有圆形标记
outer_ring_linewidth <- 0.5 # mm:外边界线宽
label_tick_length <- 0.025 # 与变量标签对齐的径向短线长度;仅用于定位,不表示数值刻度
label_tick_linewidth <- 0.45 # mm:标签定位短线的线宽
label_radius <- 1.25 # 外围变量标签的起点
font_candidates <- c("Microsoft YaHei", "PingFang SC", "Noto Sans CJK SC", "WenQuanYi Micro Hei")
palette_colors <- c(negative = "#2166AC", zero = "#F7F7F7", positive = "#B2182B")
ink <- "#283238"
required_packages <- c("ggplot2", "dplyr", "cowplot", "ragg", "systemfonts", "jsonlite", "digest")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages)) {
stop("请先安装这些包:", paste(missing_packages, collapse = ", "))
}
installed_fonts <- unique(systemfonts::system_fonts()$family)
font_family <- font_candidates[font_candidates %in% installed_fonts][1]
if (is.na(font_family)) {
stop("未找到可显示中文的字体。请把 font_candidates 改为本机已安装的中文字体。")
}
# Windows 原生绘图窗需注册字体名;此映射仅在当前 R 会话有效。
if (.Platform$OS.type == "windows") {
do.call(grDevices::windowsFonts,
setNames(list(grDevices::windowsFont(font_family)), font_family))
}
# 2. 下载官方数据,并保留原始文件 ---------------------------------------
# 每行是一套在 Ames, Iowa 于 2006—2010 年间交易的住宅;PID 为住宅编号。
# PID 按字符读入,保留前导零,且不参与相关性计算。
data_url <- "https://jse.amstat.org/v19n3/decock/AmesHousing.txt"
dictionary_url <- "https://jse.amstat.org/v19n3/decock/DataDocumentation.txt"
article_url <- "https://mp.weixin.qq.com/s/p2VFKrhRQl6nMDb22ZanPA"
data_path <- file.path(data_dir, "AmesHousing.txt")
dictionary_path <- file.path(data_dir, "DataDocumentation.txt")
if (!file.exists(data_path)) {
download.file(data_url, data_path, mode = "wb", quiet = FALSE)
}
if (!file.exists(dictionary_path)) {
download.file(dictionary_url, dictionary_path, mode = "wb", quiet = FALSE)
}
ames_raw <- read.delim(
data_path, check.names = FALSE, stringsAsFactors = FALSE,
colClasses = "character", na.strings = c("NA", ""), fileEncoding = "UTF-8"
)
stopifnot(nrow(ames_raw) == 2930L, ncol(ames_raw) == 82L)
stopifnot(!anyNA(ames_raw$PID), !anyDuplicated(ames_raw$PID))
# 变量选择服务于教学展示:售价、建成年份与同样的 16 个房屋特征进行比较。
# 未用售价筛选住宅,未按 p 值选变量,也未删除零值或较大的合法值。
targets <- data.frame(
variable = c("SalePrice", "Year Built"),
label = c("售价", "建成年份"),
unit = c("美元", "年"),
half = c("upper", "lower"),
stringsAsFactors = FALSE
)
features <- data.frame(
variable = c(
"Mas Vnr Area", "BsmtFin SF 1", "Bsmt Unf SF", "Total Bsmt SF",
"1st Flr SF", "2nd Flr SF", "Bedroom AbvGr", "TotRms AbvGrd",
"Full Bath", "Half Bath", "Garage Cars", "Garage Area",
"Wood Deck SF", "Open Porch SF", "Fireplaces", "Year Remod/Add"
),
label = c(
"饰面面积", "地下室已装面积", "地下室未装面积", "地下室总面积",
"一层面积", "二层面积", "地上卧室数", "地上房间数",
"全卫数", "半卫数", "车库车位", "车库面积",
"木平台面积", "开放门廊面积", "壁炉数", "改建年份"
),
definition = c(
"外墙砖石饰面面积", "第一类已装修地下室面积", "未装修地下室面积",
"地下室总面积", "一层面积", "二层面积", "地下室以上卧室数",
"地下室以上房间总数(不含浴室)", "地下室以上完整浴室数",
"地下室以上半浴室数", "车库容量(车位数)", "车库面积",
"木制露台面积", "开放式门廊面积", "壁炉数", "改建或增建年份"
),
unit = c(rep("平方英尺", 6), rep("个", 4), "车位", rep("平方英尺", 3), "个", "年"),
stringsAsFactors = FALSE
)
all_variables <- c(targets$variable, features$variable)
stopifnot(length(all_variables) == 18L, !anyDuplicated(all_variables))
stopifnot(all(all_variables %in% names(ames_raw)))
ames_numeric <- ames_raw[, all_variables, drop = FALSE]
for (variable in all_variables) {
converted <- as.numeric(ames_numeric[[variable]])
# 转换后若出现新的 NA,说明存在无法解析的值,应先查明,不能静默删除。
stopifnot(!any(is.na(converted) & !is.na(ames_numeric[[variable]])))
ames_numeric[[variable]] <- converted
}
complete_rows <- complete.cases(ames_numeric)
ames_complete <- ames_numeric[complete_rows, , drop = FALSE]
n_total <- nrow(ames_raw)
n_complete <- nrow(ames_complete)
n_excluded <- sum(!complete_rows)
stopifnot(n_complete == 2905L, n_excluded == 25L)
stopifnot(all(vapply(ames_complete, function(x) all(is.finite(x)), logical(1))))
stopifnot(all(vapply(ames_complete, function(x) length(unique(x)) > 1, logical(1))))
# 3. 计算 Spearman 相关及 32 项检验的 BH 校正 ---------------------------
# 房价、面积等分布偏斜;Spearman 比较秩次之间的单调关系。
# 计数和年份有大量同值,exact = FALSE 使用近似检验。
# 先对全部 32 项检验一次性校正,再加星号;校正不改变相关系数 rho。
correlation_results <- do.call(rbind, lapply(seq_len(nrow(targets)), function(i) {
do.call(rbind, lapply(seq_len(nrow(features)), function(j) {
result <- cor.test(
ames_complete[[targets$variable[i]]],
ames_complete[[features$variable[j]]],
method = "spearman", exact = FALSE, alternative = "two.sided"
)
data.frame(
target = targets$variable[i], target_label = targets$label[i],
feature = features$variable[j], feature_label = features$label[j],
feature_order = j, half = targets$half[i],
n = n_complete, rho = unname(result$estimate), p_value = result$p.value,
stringsAsFactors = FALSE
)
}))
}))
correlation_results$p_adj <- p.adjust(correlation_results$p_value, method = "BH")
correlation_results$significance <- dplyr::case_when(
correlation_results$p_adj < 0.001 ~ "***",
correlation_results$p_adj < 0.01 ~ "**",
correlation_results$p_adj < 0.05 ~ "*",
TRUE ~ ""
)
stopifnot(nrow(correlation_results) == 32L)
stopifnot(all(abs(correlation_results$rho) <= 1))
stopifnot(all(correlation_results$p_adj >= 0 & correlation_results$p_adj <= 1))
# 独立算一次秩相关矩阵,核对所有 32 个系数。
correlation_reference <- cor(ames_complete, method = "spearman")
reference_values <- mapply(
function(target, feature) correlation_reference[target, feature],
correlation_results$target, correlation_results$feature
)
stopifnot(isTRUE(all.equal(correlation_results$rho, as.numeric(reference_values))))
# 4. 构建圆周坐标与配色 -------------------------------------------------
# 一个圆的上半圈为售价,下半圈为建成年份;房屋特征上下镜像排列。
# 圆形标记内的数值为相关系数 rho,星号为 BH 校正后的 p 值等级。
# 标记固定大小;半径与圆周顺序均不编码相关强弱。
# 半步偏移使左右两端的标记分开,避免落在半圈分界线上。
upper_angles <- 180 - (seq_len(nrow(features)) - 0.5) * 180 / nrow(features)
lower_angles <- 360 - upper_angles
plot_data <- correlation_results %>%
mutate(
angle_deg = ifelse(half == "upper", upper_angles[feature_order], lower_angles[feature_order]),
angle_rad = angle_deg * pi / 180,
x = marker_center_radius * cos(angle_rad), y = marker_center_radius * sin(angle_rad),
label_x = label_radius * cos(angle_rad),
label_y = label_radius * sin(angle_rad),
rho_label = sub("^-", "−", sprintf("%.2f", rho)),
label_angle = ifelse(
angle_deg > 90 & angle_deg < 270, angle_deg - 180,
ifelse(angle_deg >= 270, angle_deg - 360, angle_deg)
),
label_hjust = ifelse(angle_deg > 90 & angle_deg < 270, 1, 0)
)
# 根据实际填色,选择对比度更高的白色或深色文字;数值和星号使用同一文字色。
fill_fun <- scales::gradient_n_pal(unname(palette_colors), values = c(0, 0.5, 1), space = "Lab")
plot_data$fill_hex <- fill_fun((plot_data$rho + 1) / 2)
relative_luminance <- function(colors) {
rgb <- t(grDevices::col2rgb(colors)) / 255
linear <- ifelse(rgb <= 0.04045, rgb / 12.92, ((rgb + 0.055) / 1.055)^2.4)
drop(linear %*% c(0.2126, 0.7152, 0.0722))
}
fill_luminance <- relative_luminance(plot_data$fill_hex)
marker_ink <- "#000000"
ink_luminance <- relative_luminance(marker_ink)
contrast_white <- 1.05 / (fill_luminance + 0.05)
contrast_ink <- (fill_luminance + 0.05) / (ink_luminance + 0.05)
plot_data$text_color <- ifelse(contrast_white > contrast_ink, "white", marker_ink)
plot_data$text_contrast <- pmax(contrast_white, contrast_ink)
stopifnot(min(plot_data$text_contrast) >= 4.5)
circle_angle <- seq(0, 2 * pi, length.out = 361)
circle_outline <- data.frame(x = outer_ring_radius * cos(circle_angle),
y = outer_ring_radius * sin(circle_angle))
# 32 条径向短线与圆内标记、外围变量标签处于相同角度,帮助逐项对应。
# 这些是标签定位线,不是相关系数刻度;它们的长度不表示统计量。
label_ticks <- plot_data %>%
mutate(
x = outer_ring_radius * cos(angle_rad),
y = outer_ring_radius * sin(angle_rad),
xend = (outer_ring_radius + label_tick_length) * cos(angle_rad),
yend = (outer_ring_radius + label_tick_length) * sin(angle_rad)
)
stopifnot(outer_ring_radius + label_tick_length < label_radius)
# 用数据坐标绘制等半径圆,图导出时圆形的尺寸与间距不会受点大小单位影响。
# 相邻圆保留空隙,外边界位于所有圆的外侧,不穿过圆形标记。
adjacent_center_distance <- 2 * marker_center_radius * sin(pi / nrow(plot_data))
stopifnot(2 * marker_radius < adjacent_center_distance)
stopifnot(marker_center_radius + marker_radius < outer_ring_radius)
marker_angles <- seq(0, 2 * pi, length.out = 81)
marker_polygons <- do.call(rbind, lapply(seq_len(nrow(plot_data)), function(i) {
data.frame(
marker_id = i,
x = plot_data$x[i] + marker_radius * cos(marker_angles),
y = plot_data$y[i] + marker_radius * sin(marker_angles),
rho = plot_data$rho[i]
)
}))
make_circle_panel <- function() {
panel_data <- plot_data
center_labels <- data.frame(
x = 0, y = c(0.27, -0.27), label = targets$label,
units = paste0("(", targets$unit, ")"), stringsAsFactors = FALSE
)
panel_plot <- ggplot() +
geom_path(data = circle_outline, aes(x, y),
colour = "#8B9597", linewidth = outer_ring_linewidth) +
geom_segment(data = label_ticks, aes(x, y, xend = xend, yend = yend),
colour = "#8B9597", linewidth = label_tick_linewidth, lineend = "round") +
annotate("segment", x = -outer_ring_radius, xend = outer_ring_radius, y = 0, yend = 0,
colour = "#CDD2D3", linewidth = 0.25, linetype = "dashed") +
geom_polygon(data = marker_polygons, aes(x, y, group = marker_id, fill = rho),
colour = "#A5ADAF", linewidth = 0.18) +
geom_text(data = panel_data,
aes(x, y = y + if (show_significance) 0.024 else 0,
label = rho_label, colour = text_color),
family = font_family, size = 3.2, fontface = "bold") +
geom_text(data = panel_data,
aes(label_x, label_y, label = feature_label, angle = label_angle, hjust = label_hjust),
family = font_family, size = 4.0, colour = ink) +
geom_text(data = center_labels, aes(x, y, label = label),
family = font_family, size = 5.4, fontface = "bold", colour = ink) +
geom_text(data = center_labels, aes(x, y = y - 0.125, label = units),
family = font_family, size = 3.8, colour = "#5D686A") +
scale_fill_gradient2(
low = palette_colors["negative"], mid = palette_colors["zero"],
high = palette_colors["positive"], midpoint = 0, limits = c(-1, 1),
breaks = c(-1, -0.5, 0, 0.5, 1), name = "Spearman ρ",
guide = guide_colourbar(
direction = "horizontal", title.position = "top", title.hjust = 0.5,
barwidth = grid::unit(7.5, "cm"), barheight = grid::unit(0.30, "cm")
)
) +
coord_fixed(xlim = c(-1.80, 1.80), ylim = c(-1.80, 1.80), expand = FALSE, clip = "off") +
theme_void(base_family = font_family) +
theme(
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(4, 4, 4, 4),
legend.position = "bottom", legend.text = element_text(size = 11, colour = ink),
legend.title = element_text(size = 12, colour = ink)
)
if (show_significance) {
panel_plot <- panel_plot +
geom_text(data = panel_data,
aes(x, y = y - 0.038, label = significance, colour = text_color),
family = font_family, size = 2.6, fontface = "bold")
}
panel_plot + scale_colour_identity()
}
circle_panel <- make_circle_panel()
# cowplot 在组装时也要测量字体;agg 支持本机中文字体,避免默认 PDF 字体警告。
# 只在本段临时切换,组装结束恢复调用者原先的设备设置。
previous_null_device <- cowplot::set_null_device("agg")
shared_legend <- cowplot::get_legend(circle_panel)
legend_plot <- cowplot::ggdraw(shared_legend)
significance_plot <- cowplot::ggdraw()
if (show_significance) {
significance_plot <- significance_plot +
cowplot::draw_label("BH 校正:* p < 0.05 ** p < 0.01 *** p < 0.001",
x = 0.5, y = 0.5,
fontfamily = font_family, size = 10.5, colour = ink)
}
circular_correlation_plot <- cowplot::plot_grid(
circle_panel + theme(legend.position = "none"), legend_plot, significance_plot,
ncol = 1, rel_heights = c(1, 0.075, 0.025)
)
cowplot::set_null_device(previous_null_device)
# 5. 导出图、实际计算结果和数据说明 -------------------------------------
ggsave(
file.path(output_dir, "ames_circular_correlation.png"), circular_correlation_plot,
device = ragg::agg_png, width = figure_width, height = figure_height,
units = "in", dpi = figure_dpi, bg = "white"
)
ggsave(
file.path(output_dir, "ames_circular_correlation.pdf"), circular_correlation_plot,
device = grDevices::cairo_pdf, width = figure_width, height = figure_height,
units = "in", bg = "white"
)
write.csv(correlation_results,
file.path(output_dir, "ames_correlation_results.csv"), row.names = FALSE, fileEncoding = "UTF-8")
variable_dictionary <- rbind(
data.frame(role = "target", variable = targets$variable, label = targets$label,
definition = targets$label, unit = targets$unit),
data.frame(role = "feature", variable = features$variable, label = features$label,
definition = features$definition, unit = features$unit)
)
write.csv(variable_dictionary,
file.path(output_dir, "ames_variable_dictionary.csv"), row.names = FALSE, fileEncoding = "UTF-8")
missing_summary <- data.frame(variable = all_variables, n_missing = colSums(is.na(ames_numeric)))
write.csv(missing_summary,
file.path(output_dir, "ames_missing_summary.csv"), row.names = FALSE, fileEncoding = "UTF-8")
write.csv(data.frame(PID = ames_raw$PID[!complete_rows]),
file.path(output_dir, "ames_excluded_pid.csv"), row.names = FALSE, fileEncoding = "UTF-8")
caption <- paste0(
"单圆环形相关性图。上半圈为售价,下半圈为建成年份;两个半圈比较同样的 16 项房屋特征。",
"等大圆形标记完整置于外边界圆环内;填色与圆内数值表示 Spearman ρ,色带范围固定为 −1 至 1;标记大小与圆周位置不表示相关强弱。",
"外环短线帮助对应圆内标记与外围变量标签,仅用于定位,不表示数值刻度。",
if (show_significance) "32 项检验一起进行 BH 校正;* p < 0.05、** p < 0.01、*** p < 0.001。" else "本图不显示显著性星号。",
"原表 2930 套住宅;所选 18 列有缺失的 25 行排除,统一使用 2905 行完整记录,零值保留。"
)
limitations <- c(
"本图用于多变量相关性展示与教学探索,不表示因果关系。",
"ρ 是未经协变量调整的 Spearman 秩相关;BH 校正仅调整 p 值,不改变 ρ。",
"Spearman 的常规 p 值检验假定观测独立;本例未对住宅的空间依赖进行调整。",
"资料覆盖 Ames 市 2006—2010 年住宅交易,不能视为全美国住宅的随机样本。",
"售价使用原表交易美元价格,未作通胀调整;这些相关也可能反映成交年份、房屋结构与地段等差异。",
"相关结果较多达到小 p 值,不代表所有关系均有较大效应;判断强弱应看 ρ。",
"本图集中展示两个核心指标与住房特征的关系,不是所选 18 个变量的完整两两相关矩阵。",
"若需要比较任意两个变量,常规矩阵热图更便于逐项读取。"
)
writeLines(enc2utf8(c(caption, "", limitations, "", paste0("数据:", data_url),
paste0("变量文档:", dictionary_url), paste0("布局参考:", article_url))),
file.path(output_dir, "caption_and_notes.txt"), useBytes = TRUE)
alt_text <- paste0(
"一个圆图,上下两半各显示 16 个固定大小的圆形标记,全部位于外边界圆环内。",
"圆心分别标注售价和建成年份;",
"外环有 32 条与变量标签对齐的径向短线,用于定位。",
"外围标签为房屋特征。蓝色表示负相关,近白色表示弱相关,红色表示正相关。",
"圆形标记内分两行标注两位小数的相关系数与显著性星号;具体 32 项系数、样本数和调整后 p 值见配套 CSV。"
)
writeLines(enc2utf8(alt_text), file.path(output_dir, "alt_text.txt"), useBytes = TRUE)
metadata <- list(
created_at = format(Sys.time(), "%Y-%m-%d %H:%M:%S %z"),
audience = "R 语言教学;多变量相关性图进阶示例",
data_source = data_url, dictionary_source = dictionary_url, layout_reference = article_url,
source_citation = "De Cock, D. (2011). Ames, Iowa: Alternative to the Boston Housing Data as an End of Semester Regression Project. Journal of Statistics Education, 19(3).",
observation_unit = "住宅;PID 唯一,2006—2010 年 Ames 市住宅交易资料",
source_rows = n_total, source_columns = ncol(ames_raw), complete_rows = n_complete,
excluded_missing_rows = n_excluded, variables_used = all_variables,
zero_values_retained = TRUE, outlier_filtering = "none",
method = "Spearman rho; two-sided cor.test; exact=FALSE",
multiple_testing = "BH adjustment across all 32 tests",
n_tests = nrow(correlation_results),
n_adjusted_p_below_005 = sum(correlation_results$p_adj < 0.05),
show_significance = show_significance,
significance_thresholds = list(`*` = 0.05, `**` = 0.01, `***` = 0.001),
significance_counts = list(
not_significant = sum(correlation_results$significance == ""),
`*` = sum(correlation_results$significance == "*"),
`**` = sum(correlation_results$significance == "**"),
`***` = sum(correlation_results$significance == "***")
),
font_family = font_family, palette_hex = as.list(palette_colors),
color_limits = c(-1, 1), layout = "single circle; upper SalePrice, lower Year Built",
marker = "fixed-radius circle; rho and stars on separate lines; fully inside outer ring",
marker_center_radius = marker_center_radius, marker_radius = marker_radius,
outer_ring_radius = outer_ring_radius,
outer_ring_linewidth_mm = outer_ring_linewidth,
label_tick_length = label_tick_length,
label_tick_linewidth_mm = label_tick_linewidth,
label_tick_count = nrow(label_ticks),
label_ticks_semantics = "label alignment guides; not a numerical scale",
label_radius = label_radius,
adjacent_marker_gap = adjacent_center_distance - 2 * marker_radius,
marker_to_ring_gap = outer_ring_radius - marker_center_radius - marker_radius,
minimum_text_contrast = min(plot_data$text_contrast),
requested_inches = c(width = figure_width, height = figure_height),
png_pixels = c(width = figure_width * figure_dpi, height = figure_height * figure_dpi),
dpi = figure_dpi, source_sha256 = digest::digest(data_path, algo = "sha256", file = TRUE),
limitations = limitations,
package_versions = setNames(lapply(required_packages, function(x) as.character(packageVersion(x))),
required_packages)
)
jsonlite::write_json(metadata, file.path(output_dir, "metadata.json"),
pretty = TRUE, auto_unbox = TRUE, null = "null")
writeLines(capture.output(sessionInfo()), file.path(output_dir, "session_info.txt"))
message("数据:", n_total, " 行;完整记录:", n_complete,
";缺失剔除:", n_excluded, ";相关比较:", nrow(correlation_results))
message("BH 校正 p < 0.05:", sum(correlation_results$p_adj < 0.05), "/32;输出:", output_dir)
# RStudio 运行后可直接查看;也可继续修改 show_significance 和 palette_colors。
if (interactive()) print(circular_correlation_plot)解释边界。 图中系数未控制成交年份、地段等因素,检验也未调整住宅之间的空间依赖,因此不能据此判断因果关系。本图展示两个核心指标与住房特征的 32 组关系;需要逐项比较所有变量时,完整相关矩阵更便于查阅。
布局思路参考公众号「农学日记」:原文链接。本例更换数据,并重新计算相关系数与检验结果。
图形应交代观测对象、变量和统计口径。读者看坐标、图例与图注,就应知道图中比较了什么、哪些差异可以直接比较。本节按选图、尺度、配色、布局、标注和出图检查的顺序展开;地图投影与空间尺度留待后续课程讨论。
同一组数据可以用于比较水平、观察变化或展示分布。先看每条记录包含哪些变量,再明确图形要回答的问题。
| 数据与问题 | 可先考虑 | 读图时核对 |
|---|---|---|
| 一列数值:看分布 | 直方图或核密度图 | 直方图的组距与纵轴口径;密度图的平滑程度 |
| 两列数值:看关系 | 散点图;点密集时可用二维分箱 | 一个点代表什么;点是否遮挡,分箱的颜色表示什么 |
| 数值随时间变化 | 折线图;非负数量需强调数值大小时可用面积图 | 相邻点是否属于同一序列;时间间隔与缺测是否清楚 |
| 数值 × 类别:比较组别 | 同一分析条件下,每类有多个对象的原始值时用箱线图、小提琴图或散点;每类仅有一个汇总值时用柱状图或点图 | 每类样本量;汇总值是总量、均值还是比例 |
| 类别 × 类别:看交叉频数或构成 | 分组柱状图、堆叠柱状图或比例热力图 | 单元格表示人数还是比例;比例以哪一组为分母 |
类别 × 类别的交叉频数,可由每条记录中的两列类别标签计数得到,原始数据无须另有一列数值。
网络图要有明确的连接。 若两列分别记录关系的起点和终点,它们可以构成边表;若两列只是同一对象的两个类别属性,则不能据此画出对象之间的关系。Q–Q 图用于比较样本与正态分布等参照分布,是分布诊断,不等同于统计检验。
直方图应交代组距,说明纵轴是频数、比例还是密度;二维分箱还需交代两个轴的格宽,以及填色表示计数还是密度。时间序列存在缺测或中断时,应断开线段或明确标注,避免连线暗示中间时点已有观测。
局部放大与限制数据范围。
coord_cartesian() 只改变视野;xlim() 或
ylim()
可能移除范围外的观测,进而改变拟合线等统计图层。使用截取范围时,应让坐标刻度清楚显示这一点。
颜色应对应清楚的变量,并由图例或直接标注说明。色板的具体设置见 5.0。
灰度打印或投影时,可补充点形、线型或纹理;辅助编码应与图例一致。
连续数值、相对基准的差值与无序类别的配色示意
facet_wrap() 或
facet_grid()。若要比较绝对水平,仍需核对各面板的坐标范围。组别分面、重叠散点处理与堆叠柱状图的替代布局;六边形分箱仅为示意,不能据此读取具体数值
轴标签应说明变量含义与单位;数据文档没有说明的口径,不要自行补写。发生对数变换、标准化或指数化时,还须交代处理方式。本章的
gdpPercap 可标为“人均 GDP”。若使用
scale_x_log10(),可写成“人均
GDP(对数刻度)”:刻度仍显示原始数值,轴上的距离对应数值的倍数变化。
代码变量名不宜原样放在图上。例如,life_exp_w
可标为“人口加权平均预期寿命(岁)”;分组缩写应在图例或图注中解释。多图应使用一致的字族和同类元素样式,轴标题、刻度和图例仍需有清楚的字号层级。
图形应先保证数据表达准确、读图判断清楚,再调整配色、字体与留白。同一组图保持一致的视觉层级,便于读者比较结果。
(本节以非地图类图形为例。)
图形完成后,缩到网页或论文中的实际尺寸,检查图例、标签和重点是否仍能看清;同一变量在不同图中尽量沿用相同颜色。
区分类别与重点:类别颜色应彼此可辨;强调色只用于需要读者注意的部分。大面积高饱和填色容易掩盖文字和细节。
表达数值结构:连续数值需要有序色带,基准两侧的偏离可用发散色带;图例应说明颜色与数值的对应关系。
保持跨图一致:同一变量在一组图中沿用相同颜色,读者无需反复辨认图例。
顺序色带:用于连续变量。
颜色应随数值呈现清楚的明度顺序,并在图例中标明数值方向。
(不建议使用“彩虹色带”,其亮度跳跃容易制造伪边界,干扰梯度判断)
发散色带:用于有明确基准值的数据(如相对某一标准的差值)。
基准值放在色带中点,两侧用不同色相表示低于和高于基准;图例应标出中点的实际数值。
(中点需有分析意义,不能仅因数据范围居中就设为基准。)
离散色板:用于分类变量。
类别之间需可区分,同时避免某一颜色过于“抢眼”导致视觉权重失衡。
(当存在“大组 + 子组”结构,可采用“组间不同色相 +
组内同色系明暗”以提示隶属关系)
下面各卡展示两套代表性色盘,便于横向比较。更多色盘和原有示例代码见下方折叠区。
按数据类型选预设色盘。
连续数值的颜色有清楚的明暗层次。
提供多套适合类别区分的期刊风格色组。
汇集色盘,并提供色觉可辨性检查。
提供 CARTO 的顺序、发散和分类色盘。
用 HCL 色彩空间生成和调整色阶。
package_palette_examples <- list(
brewer = list(
"Set2 · 分类" = RColorBrewer::brewer.pal(7, "Set2"),
"RdBu · 发散" = RColorBrewer::brewer.pal(7, "RdBu")
),
viridis = list(
"Viridis · 连续" = viridisLite::viridis(7),
"Cividis · 连续" = viridisLite::cividis(7)
),
ggsci = list(
"NPG · 分类" = ggsci::pal_npg()(7),
"Lancet · 分类" = ggsci::pal_lancet()(7)
),
cols4all = list(
"Tableau Color Blind · 分类" = cols4all::c4a("tableau.color_blind", n = 7),
"Tableau Orange Gold · 连续" = cols4all::c4a("tableau.orange_gold", n = 7)
),
carto = list(
"Burg · 连续" = rcartocolor::carto_pal(7, "Burg"),
"Mint · 连续" = rcartocolor::carto_pal(7, "Mint")
),
colorspace = list(
"Blues 3 · 连续" = colorspace::sequential_hcl(7, palette = "Blues 3"),
"Blue-Red 3 · 发散" = colorspace::diverging_hcl(7, palette = "Blue-Red 3")
)
)RColorBrewer
推荐访问viridis
推荐访问viridis色盘与名称
ggsci
推荐访问# install.packages('ggsci')
library("ggsci")
library("ggplot2")
library("gridExtra")
p1 <- example_scatterplot()
p2 <- example_barplot()
p1_npg <- p1 + scale_color_npg()
p2_npg <- p2 + scale_fill_npg()
grid.arrange(p1_npg, p2_npg, ncol = 2) # 展示 NPG 系列配色cols4all
推荐访问# install.packages('cols4all')
library(cols4all)
# cols4all 汇集多种色盘,并提供色觉差异相关的检查工具。
# 可在本地 RStudio 运行 c4a_gui() 浏览色盘。
# c4a_gui()
# 这里使用 tableau.orange_gold;可用 c4a_palettes() 查询其他色盘名称。
pal_colors <- c4a("tableau.orange_gold", n = 10)
c4a_plot(pal_colors) # 预览该色盘rcartocolor
推荐访问colorspace
推荐访问下面三个网站分别适合试配色、比较图形效果和按数据类型筛选色盘。截图可放大查看。
其他调色工具 Coolors(快速组合)、Adobe Color(色轮)、Paletton(配色关系)、Colormind(自动生成)。
取色与灵感 Flat UI Color Picker、Brand Colors、Color Leap。
延伸阅读 Graphics and Colors、R 图形教程。
下列影视画面附有提取出的色盘,可对照观察主色、辅色和强调色的搭配及使用比例。
《小丑》
《安娜·卡列尼娜》
《楚门的世界》
《布达佩斯大饭店》
《鱿鱼游戏》
《美国人质》
《朱诺》
《西西里的美丽传说》
如何借鉴这些色盘
影视作品的颜色主要服务于画面表达,适合用来学习搭配。学术绘图还要让读者准确辨认数据:直接照搬色盘,可能出现类别区分不清、文字对比不足,或在色觉差异、灰度打印等条件下难以辨认的问题。
可先从作品中选两三种颜色,再结合数据含义调整。连续数值需要色阶时,可以用渐变工具补足中间色,并检查明暗顺序;类别数据则要逐一挑选容易区分的颜色,不宜简单插值,必要时还可用点形或线型辅助。最后将颜色放入实际图表,检查标签、图例和灰度效果。
色值可以手动输入,也可以从配色包取得。写进图表时,先看颜色是整层固定,还是随类别或连续数值变化。下面三图均使用 Gapminder 样本数据。
HEX 色值写作 #RRGGBB。若需透明度,可用
scales::alpha("#4C78A8", 0.7);第二个参数表示不透明度,范围为
0–1。
geom_bar(fill = "#4C78A8")
在 aes() 外写入色值,所有柱形共用这一种颜色。
aes(fill = continent)scale_fill_manual(values =
my_cols)
用类别名固定颜色。色值可手写,也可从 RColorBrewer 或
ggsci 取出。
aes(fill = avg_life)scale_fill_viridis_c()
色阶表示各大洲样本国家预期寿命的等权均值。
library(ggplot2)
library(dplyr)
library(gapminder)
df_2007 <- gapminder %>% filter(year == 2007)
continent_zh <- c(
Africa = "非洲", Americas = "美洲", Asia = "亚洲",
Europe = "欧洲", Oceania = "大洋洲"
)
ggplot(df_2007, aes(x = continent)) +
geom_bar(fill = "#4C78A8", width = 0.68) +
coord_flip() +
scale_x_discrete(labels = continent_zh) +
labs(x = NULL, y = "样本国家数") +
theme_minimal(base_size = 10) +
theme(panel.grid.major.y = element_blank())library(ggplot2)
library(dplyr)
library(gapminder)
df_2007 <- gapminder %>% filter(year == 2007)
continent_zh <- c(
Africa = "非洲", Americas = "美洲", Asia = "亚洲",
Europe = "欧洲", Oceania = "大洋洲"
)
my_cols <- c(
Africa = "#4C78A8", Americas = "#F58518",
Asia = "#54A24B", Europe = "#E45756", Oceania = "#72B7B2"
)
# 也可从包内取色,保留类别名与颜色的对应关系:
# my_cols <- setNames(RColorBrewer::brewer.pal(5, "Set2"), levels(df_2007$continent))
# my_cols <- setNames(ggsci::pal_npg()(5), levels(df_2007$continent))
ggplot(df_2007, aes(x = continent, fill = continent)) +
geom_bar(width = 0.68) +
coord_flip() +
scale_x_discrete(labels = continent_zh) +
scale_fill_manual(values = my_cols, guide = "none") +
labs(x = NULL, y = "样本国家数") +
theme_minimal(base_size = 10) +
theme(panel.grid.major.y = element_blank())library(ggplot2)
library(dplyr)
library(gapminder)
continent_zh <- c(
Africa = "非洲", Americas = "美洲", Asia = "亚洲",
Europe = "欧洲", Oceania = "大洋洲"
)
life_by_year <- gapminder %>%
group_by(continent, year) %>%
summarise(avg_life = mean(lifeExp), .groups = "drop") # 样本国家等权平均
ggplot(life_by_year, aes(x = year, y = continent, fill = avg_life)) +
geom_tile(width = 4.9, height = 0.9) +
scale_x_continuous(breaks = c(1952, 1972, 1992, 2007)) +
scale_y_discrete(labels = continent_zh) +
scale_fill_viridis_c(name = "预期寿命\n均值(岁)") +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 9) +
theme(panel.grid = element_blank(), legend.key.height = grid::unit(0.7, "cm"))若用 viridis 给类别着色,可改用
scale_fill_viridis_d(),并检查类别之间是否容易区分。
以下六套配色使用相同的数据与六图模板,便于比较颜色在不同图形中的表现。后三套色值据课堂提供的参考图整理,所有图形均由本章 R 代码生成。选色时仍需检查数据含义和可读性;点击图片可放大查看细节与色号。
图 a、b、e、f 使用固定随机种子的模拟数据;图 c 使用
ggplot2::diamonds;图 d 展示颜色过渡。
紫青黄渐变(7 色)
紫青黄渐变在六种图形上的配色预览
柔和对比(5 色)
柔和对比在六种图形上的配色预览
暮紫暖金(5 色)
暮紫暖金在六种图形上的配色预览
黄昏绮景(7 色)
黄昏绮景在六种图形上的配色预览
自然活力(6 色)
自然活力在六种图形上的配色预览
经典碰撞(5 色)
经典碰撞在六种图形上的配色预览
要单独试色,可在脚本开头的 palette_sets 中添加 HEX
色值,将 selected_palette 改为新方案的
ID,再运行脚本末尾的两行示例代码。本页固定展示六套方案,便于直接对照。
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(cowplot)
})
# 色盘 ID 统一为“颜色特征_颜色数”;单独运行脚本时修改 selected_palette 即可试色。
palette_sets <- list(
purple_teal_yellow_7 = c(
"#44035b", "#404185", "#31688e", "#1f918d",
"#38b775", "#90d543", "#f8e620"
),
soft_contrast_5 = c(
"#fbf1c2", "#566970", "#e67a7f", "#9d3151", "#87baa7"
),
mauve_gold_5 = c(
"#3E2564", "#7B5E8E", "#B898A0", "#D4B896", "#E5C87B"
),
dusk_teal_coral_7 = c(
"#274753", "#297270", "#299D8F", "#8AB07C",
"#E7C66B", "#F3A361", "#E66D50"
),
vivid_nature_6 = c(
"#A30543", "#F36F43", "#FBDA83", "#E9F4A3",
"#80CBA4", "#4965B0"
),
red_blue_contrast_5 = c(
"#7A0101", "#BE1420", "#FBF0D9", "#012F48", "#669ABA"
),
# 两套已有试配色保留在脚本中,供后续替换展示方案。
navy_mauve_amber_5 = c(
"#1E3E6E", "#5A6E9E", "#9E8EAA", "#D49E8A", "#F0B878"
),
teal_lilac_5 = c(
"#1A3A3A", "#2E6E6E", "#5A8E9E", "#8E6E9E", "#B888C4"
)
)
selected_palette <- "purple_teal_yellow_7"
# 传入 5 色及以上 HEX 向量,返回可直接 print() 或 ggsave() 的 2 行 3 列六图版式。
make_palette_board <- function(palette_colors) {
if (length(palette_colors) < 5L || anyNA(palette_colors)) {
stop("palette_colors 至少需要 5 个有效颜色。")
}
palette_colors <- as.character(palette_colors)
invisible(grDevices::col2rgb(palette_colors))
pal <- palette_colors
pick_palette <- function(colors, n) {
colors[round(seq(1, length(colors), length.out = n))]
}
ink <- "#283238"
canvas_bg <- "white"
panel_theme <- theme_bw(base_size = 11, base_family = "sans") +
theme(
panel.grid = element_blank(),
panel.border = element_rect(colour = "#4d5556", linewidth = 0.35),
axis.line = element_blank(),
axis.ticks = element_line(colour = ink, linewidth = 0.35),
axis.text = element_text(colour = ink, size = 9),
axis.title = element_text(colour = ink, face = "bold", size = 11),
legend.text = element_text(colour = ink, size = 8),
legend.key.size = grid::unit(0.33, "cm"),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(8, 9, 8, 9)
)
# 第一行左侧:固定随机种子的模拟簇,用于展示多类别散点
set.seed(20261002)
centers <- data.frame(
group = LETTERS[1:8],
cx = c(5, 5, 9, 9, 13, 14, 18, 19),
cy = c(5, 10, 4, 8, 7, 12, 6, 9)
)
scatter_df <- do.call(rbind, lapply(seq_len(nrow(centers)), function(i) {
data.frame(
group = centers$group[i],
x = rnorm(130, centers$cx[i], if (i %% 2 == 0) 1.35 else 1.05),
y = rnorm(130, centers$cy[i], 1.15)
)
}))
cluster_cols <- setNames(
if (length(pal) == 5L) pal[c(1, 2, 2, 3, 3, 4, 4, 5)]
else grDevices::colorRampPalette(pal)(8),
LETTERS[1:8]
)
p_scatter <- ggplot(scatter_df, aes(x, y, fill = group)) +
geom_point(shape = 21, colour = ink, stroke = 0.12,
size = 1.15, alpha = 0.82) +
geom_label(
data = centers, aes(x = cx, y = cy, label = group),
inherit.aes = FALSE, colour = ink, fill = scales::alpha("white", 0.82),
fontface = "bold", size = 3.6,
label.padding = grid::unit(0.07, "lines"), linewidth = 0
) +
scale_fill_manual(values = cluster_cols, guide = "none") +
scale_x_continuous(breaks = c(5, 10, 15, 20)) +
scale_y_continuous(breaks = c(4, 8, 12, 16)) +
coord_cartesian(xlim = c(0, 23), ylim = c(1, 16)) +
labs(x = "x", y = "y") +
panel_theme
# 第一行中间:模拟四组数值,只比较这套颜色在小提琴和箱线图上的效果
samples <- paste("Sample", 1:4)
violin_df <- data.frame(
sample = factor(rep(samples, each = 180), levels = samples),
value = unlist(Map(
function(mu, sd) rnorm(180, mu, sd),
c(54, 59, 61, 49), c(8, 9, 10, 7)
))
)
p_violin <- ggplot(violin_df, aes(sample, value, fill = sample)) +
geom_violin(trim = FALSE, width = 0.82, alpha = 0.82, colour = NA) +
geom_boxplot(
width = 0.16, fill = "white", colour = ink, linewidth = 0.42,
outlier.shape = 16, outlier.size = 0.8
) +
scale_fill_manual(values = setNames(pick_palette(pal, 4), samples), guide = "none") +
coord_cartesian(ylim = c(15, 100)) +
labs(x = "Sample", y = "Value") +
panel_theme
# 花瓣用于展示色系过渡,不对应统计量。
make_petal_panel <- function(pal) {
n <- 15L
petal_cols <- grDevices::colorRampPalette(pal)(n)
theta <- pi / 2 + (seq_len(n) - 1) * 2 * pi / n
t <- seq(0, 2 * pi, length.out = 121)
petals <- do.call(rbind, lapply(seq_len(n), function(i) {
a <- theta[i]
radial <- 1.13 + 1.12 * cos(t)
tangential <- 0.43 * sin(t)
data.frame(
id = factor(i, levels = seq_len(n)),
x = radial * cos(a) - tangential * sin(a),
y = radial * sin(a) + tangential * cos(a)
)
}))
angle_deg <- (theta * 180 / pi) %% 360
angle_deg <- ifelse(angle_deg > 90 & angle_deg < 270, angle_deg - 180, angle_deg)
labels <- data.frame(
id = seq_len(n),
x_num = 1.68 * cos(theta), y_num = 1.68 * sin(theta),
x_outer = 2.57 * cos(theta), y_outer = 2.57 * sin(theta),
angle = angle_deg
)
core <- data.frame(x = 0.72 * cos(t), y = 0.72 * sin(t))
ggplot() +
geom_polygon(data = petals, aes(x, y, group = id, fill = id),
alpha = 0.62, colour = NA) +
geom_polygon(data = core, aes(x, y), fill = "white", colour = "white") +
geom_text(data = labels, aes(x_num, y_num, label = id), size = 3.1, colour = ink) +
geom_text(data = labels,
aes(x_outer, y_outer, label = paste0("A", id), angle = angle),
size = 3.5, colour = ink) +
annotate("text", x = 0, y = 0, label = paste(length(pal), "COLORS"),
fontface = "bold", size = 3.0, colour = ink) +
scale_fill_manual(values = setNames(petal_cols, seq_len(n)), guide = "none") +
coord_equal(xlim = c(-2.95, 2.95), ylim = c(-2.95, 2.95),
expand = FALSE, clip = "off") +
theme_void(base_family = "sans") +
theme(
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(8, 8, 8, 8)
)
}
p_petals <- make_petal_panel(pal)
# 第一行右侧:真实 diamonds 数据;每个切工内显示价格五分位的占比
price_breaks <- quantile(ggplot2::diamonds[["price"]], probs = seq(0, 1, 0.2))
diamonds_df <- ggplot2::diamonds %>%
mutate(price_band = cut(
price, breaks = price_breaks, include.lowest = TRUE,
labels = paste0("Q", 1:5)
))
band_cols <- setNames(rev(pick_palette(pal, 5)), paste0("Q", 1:5))
p_stacked <- ggplot(diamonds_df, aes(cut, fill = price_band)) +
geom_bar(position = position_fill(reverse = TRUE), width = 0.66,
colour = "white", linewidth = 0.28) +
scale_fill_manual(values = band_cols, name = "Price quintile",
guide = guide_legend(nrow = 1, title.position = "top")) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1),
expand = expansion(mult = c(0, 0))) +
labs(x = "Cut", y = "Share") +
panel_theme +
theme(
axis.text.x = element_text(size = 8.5),
legend.position = "bottom",
legend.justification = "left",
legend.title = element_text(size = 8.5, face = "bold"),
legend.margin = margin(0, 0, 0, 0)
)
# 第二行中间:模拟起伏曲面;用 lattice 绘制真正的三维表面
surface_df <- expand.grid(
x = seq(-3, 3, length.out = 45),
y = seq(-3, 3, length.out = 45)
)
surface_df$z <- with(
surface_df,
sin(1.4 * x) * cos(1.1 * y) + 0.25 * cos(2 * x + y)
)
surface_path <- if (length(pal) == 5L) pal[c(2, 5, 1, 3, 4)] else pal
surface_cols <- grDevices::colorRampPalette(surface_path)(100)
surface_trellis <- lattice::wireframe(
z ~ x * y, data = surface_df,
drape = TRUE, colorkey = FALSE,
col = "#ffffff66", lwd = 0.2,
col.regions = surface_cols,
screen = list(z = 38, x = -58, y = 0),
xlab = "x", ylab = "y", zlab = "z",
scales = list(draw = FALSE),
par.settings = list(axis.line = list(col = "#98a2a3", lwd = 0.45)),
aspect = c(1, 1, 0.6)
)
surface_grob <- grid::grid.grabExpr(print(surface_trellis))
p_surface <- cowplot::ggdraw() +
theme(plot.background = element_rect(fill = "white", colour = NA)) +
cowplot::draw_grob(surface_grob, x = -0.18, y = -0.10,
width = 1.36, height = 1.22)
# 第二行右侧:模拟风向与风速,检验配色在极坐标中的效果
set.seed(20261003)
wind_n <- 900L
wind_centers <- sample(
c(30, 120, 225, 310), wind_n, replace = TRUE,
prob = c(0.18, 0.12, 0.48, 0.22)
)
wind_directions <- c(
"N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE",
"S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"
)
speed_levels <- c("0-2", "2-4", "4-6", "6-8", ">8")
wind_df <- data.frame(
direction_deg = rnorm(wind_n, wind_centers, 22) %% 360,
speed = rgamma(wind_n, shape = 2.8, scale = 1.7)
) %>%
mutate(
sector = factor(
wind_directions[floor(((direction_deg + 11.25) %% 360) / 22.5) + 1L],
levels = wind_directions
),
speed_band = cut(
speed, breaks = c(-Inf, 2, 4, 6, 8, Inf),
labels = speed_levels
)
)
rose_df <- wind_df %>%
count(sector, speed_band, .drop = FALSE) %>%
mutate(percent = 100 * n / sum(n))
max_sector_percent <- max(tapply(rose_df$percent, rose_df$sector, sum))
upper_percent <- 5 * ceiling(max_sector_percent / 5)
p_wind <- ggplot(rose_df, aes(sector, percent, fill = speed_band)) +
geom_col(
position = position_stack(reverse = TRUE),
width = 0.94, colour = "white", linewidth = 0.18
) +
coord_polar(start = -pi / 16) +
scale_x_discrete(
labels = function(x) ifelse(
x %in% c("N", "NE", "E", "SE", "S", "SW", "W", "NW"), x, ""
)
) +
scale_y_continuous(
limits = c(0, upper_percent),
breaks = seq(0, upper_percent, by = 5),
labels = function(x) paste0(x, "%"),
expand = expansion(mult = c(0, 0))
) +
scale_fill_manual(
values = setNames(pick_palette(pal, 5), speed_levels),
name = "Wind speed (sim., m/s)",
guide = guide_legend(nrow = 1, title.position = "top")
) +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 11, base_family = "sans") +
theme(
panel.grid.major.x = element_line(colour = "#e7e9e9", linewidth = 0.25),
panel.grid.major.y = element_line(colour = "#d7dcdc", linewidth = 0.35),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = ink, size = 9, face = "bold"),
axis.text.y = element_text(colour = ink, size = 7),
axis.ticks = element_blank(),
legend.position = "bottom",
legend.title = element_text(size = 8.5, face = "bold"),
legend.text = element_text(size = 8),
legend.key.size = grid::unit(0.33, "cm"),
legend.margin = margin(0, 0, 0, 0),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(8, 9, 8, 9)
)
# 下方色卡:用固定物理尺寸绘制正方形,完整显示传入色盘的 HEX 与 RGB。
rgb_text <- apply(t(grDevices::col2rgb(pal)), 1, paste, collapse = ", ")
swatches <- cowplot::ggdraw() +
theme(plot.background = element_rect(fill = "white", colour = NA))
swatch_n <- length(pal)
for (i in seq_along(pal)) {
center <- (i - 0.5) / swatch_n
swatches <- swatches +
cowplot::draw_grob(
grid::roundrectGrob(
width = grid::unit(0.44, "in"),
height = grid::unit(0.44, "in"),
r = grid::unit(0.05, "in"),
gp = grid::gpar(fill = pal[i], col = NA)
),
x = center - 0.45 / swatch_n, y = 0.39,
width = 0.9 / swatch_n, height = 0.53
) +
cowplot::draw_text(tolower(pal[i]), x = center, y = 0.31,
size = if (swatch_n <= 5L) 10 else 9) +
cowplot::draw_text(rgb_text[i], x = center, y = 0.13,
size = if (swatch_n <= 5L) 8.5 else 7.5)
}
# 编号占用独立的上沿留白,不压在图内的坐标轴或图例上。
labeled_panel <- function(plot, panel_label) {
cowplot::ggdraw() +
theme(plot.background = element_rect(fill = "white", colour = NA)) +
cowplot::draw_plot(plot, x = 0.01, y = 0.005, width = 0.98, height = 0.895) +
cowplot::draw_label(
panel_label, x = 0.035, y = 0.975,
hjust = 0, vjust = 1, size = 11,
fontface = "bold", colour = ink
)
}
# 两行三列六图 + 一行色卡
plots <- cowplot::plot_grid(
labeled_panel(p_scatter, "a)"),
labeled_panel(p_violin, "b)"),
labeled_panel(p_stacked, "c)"),
labeled_panel(p_petals, "d)"),
labeled_panel(p_surface, "e)"),
labeled_panel(p_wind, "f)"),
ncol = 3, rel_heights = c(1, 1)
)
palette_board <- cowplot::ggdraw() +
theme(plot.background = element_rect(fill = canvas_bg, colour = NA)) +
cowplot::draw_grob(
grid::roundrectGrob(
r = grid::unit(0.025, "snpc"),
gp = grid::gpar(fill = "white", col = "#e5e8e8", lwd = 0.7)
),
x = 0.025, y = 0.205, width = 0.95, height = 0.775
) +
cowplot::draw_plot(plots, x = 0.045, y = 0.215, width = 0.91, height = 0.75) +
cowplot::draw_grob(
grid::roundrectGrob(
r = grid::unit(0.025, "snpc"),
gp = grid::gpar(fill = "white", col = "#e5e8e8", lwd = 0.7)
),
x = 0.025, y = 0.02, width = 0.95, height = 0.16
) +
cowplot::draw_plot(swatches, x = 0.045, y = 0.03, width = 0.91, height = 0.14)
palette_board
}
# 使用示例(取消注释后执行):
# palette_board <- make_palette_board(palette_sets[[selected_palette]])
# print(palette_board)
本节分为三个环节:先运行两张示例图,检查数据筛选与图层设置;再用两份数据完成小组作图报告;最后通过七道题回顾本章的出图方法与读图判断。
练习使用 fivethirtyeight::hate_crimes。数据共有
51 条州级记录,包括美国 50 个州和
DC。两张图分别使用不同的变量,因此只排除各自所需列中的缺失值:散点图保留
50 条,分布图保留 51 条。不要对全部 13
列直接运行 drop_na()。
建议新建
chapter_3_visual_practice.Rmd。在正文写明一行记录代表什么,再分别插入
R 代码块运行以下两段代码。完整代码在两图后方的折叠区,可直接复制到自己的
Rmd。运行前,如本机尚未安装,可先在 Console 中一次性安装
ggplot2、dplyr、fivethirtyeight
和 scales;任务 2 还需 gapminder。
2015 年州级基尼系数与 FBI 在 2010—2015 年记录的每 10 万人年均仇恨犯罪数。两条拟合线分别使用含 DC 和不含 DC 的观测;图例给出各自的 Pearson 相关系数。
图中每个点对应一个州或 DC。只标注 DC,可以看清其位置而不让 50 个缩写互相遮挡。两条线及其相关系数显示,DC 对这组数据的描述性关系影响很大;这不能证明基尼系数导致仇恨犯罪报告数变化。
2016 年各州及 DC 的家庭收入中位数,按特朗普得票比例是否超过 50% 分组。小提琴轮廓表示平滑后的州级分布,箱体表示中位数和四分位距,深色点是全部 51 条原始记录。
小提琴轮廓来自平滑密度估计,尾端可能超出观测值范围;51 个深色点才是实际记录。箱体概括中位数和四分位距。图中比较的是各州家庭收入中位数的分布,并非各州所有家庭收入的合并分布。
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
})
# 优先使用本机已安装的中文字体;没有这些字体时回退到 sans。
# 如果中文字形异常,可把本机字体名称加入 font_candidates。
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
font_candidates <- c("Microsoft YaHei", "PingFang SC", "Noto Sans CJK SC")
installed_font <- font_candidates[
font_candidates %in% systemfonts::system_fonts()$family
]
if (length(installed_font) > 0L) font_family <- installed_font[1]
}
# 只检查本图横纵轴所需的两列,避免因其他列缺失而删去观测。
scatter_data <- fivethirtyeight::hate_crimes %>%
filter(
!is.na(gini_index),
!is.na(avg_hatecrimes_per_100k_fbi)
)
# DC 仍在 scatter_data 中;单独取出只是为了高亮和标注。
dc_data <- scatter_data %>% filter(state_abbrev == "DC")
# 分别计算含/不含 DC 的 Pearson r,检查这个点对结果的影响。
r_with_dc <- cor(
scatter_data$gini_index,
scatter_data$avg_hatecrimes_per_100k_fbi,
method = "pearson"
)
without_dc <- scatter_data %>% filter(state_abbrev != "DC")
r_without_dc <- cor(
without_dc$gini_index,
without_dc$avg_hatecrimes_per_100k_fbi,
method = "pearson"
)
# 为两条拟合线准备数据:同一条 DC 记录只在“含 DC”组出现。
# 这里比较的是对单个观测的敏感性,并不是两组独立样本的比较。
fit_data <- bind_rows(
scatter_data %>% mutate(fit_group = "含 DC"),
without_dc %>% mutate(fit_group = "不含 DC")
) %>%
mutate(fit_group = factor(fit_group, levels = c("不含 DC", "含 DC")))
# 图例中的 r 来自上方两次 cor();n 是参与各自计算的记录数。
fit_labels <- c(
sprintf("不含 DC:r = %.2f,n = %d", r_without_dc, nrow(without_dc)),
sprintf("含 DC:r = %.2f,n = %d", r_with_dc, nrow(scatter_data))
)
ggplot(
scatter_data,
aes(x = gini_index, y = avg_hatecrimes_per_100k_fbi)
) +
# 两条线均为线性拟合;se = FALSE 不绘制置信带。
# 不含 DC 的线只画在该组实际观测到的横轴范围内。
geom_smooth(
data = fit_data, aes(colour = fit_group),
method = "lm", formula = y ~ x, se = FALSE,
linewidth = 0.85
) +
geom_point(colour = "#587684", size = 2.2, alpha = 0.82) +
# 在原有点上重画 DC;仍计作一个观测,并没有增加样本量。
geom_point(data = dc_data, colour = "#A94F5B", size = 3.4) +
geom_text(
data = dc_data, aes(label = state_abbrev),
nudge_x = -0.003, nudge_y = 0.25, hjust = 1,
colour = "#913F4A", fontface = "bold", size = 3.6
) +
scale_colour_manual(
values = c("不含 DC" = "#23817D", "含 DC" = "#A94F5B"),
breaks = c("不含 DC", "含 DC"), labels = fit_labels, name = NULL
) +
# 纵轴从 0 起;上方留白让 DC 的标签不贴边。
scale_x_continuous(expand = expansion(mult = c(0.04, 0.08))) +
scale_y_continuous(
limits = c(0, NA),
breaks = seq(0, 12, by = 3),
expand = expansion(mult = c(0, 0.08))
) +
labs(
x = "基尼系数(2015 年)",
y = "FBI 年均报告数(每 10 万人)"
) +
theme_classic(base_size = 11, base_family = font_family) +
theme(
panel.grid.major.y = element_line(colour = "#E9EDEF", linewidth = 0.3),
axis.line = element_line(colour = "#6F7B80", linewidth = 0.3),
axis.ticks = element_line(colour = "#6F7B80", linewidth = 0.3),
axis.text = element_text(colour = "#33464E"),
axis.title = element_text(colour = "#243941"),
legend.position = "top",
legend.justification = "left",
legend.text = element_text(size = 9, colour = "#33464E"),
legend.key.width = grid::unit(1.7, "lines"),
plot.margin = margin(8, 12, 8, 8)
)suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
})
# 每段代码都能单独复制运行,因此在这里重复设置字体。
# 如本机没有列出的字体,可改 font_candidates,或让 sans 自动回退。
font_family <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
font_candidates <- c("Microsoft YaHei", "PingFang SC", "Noto Sans CJK SC")
installed_font <- font_candidates[
font_candidates %in% systemfonts::system_fonts()$family
]
if (length(installed_font) > 0L) font_family <- installed_font[1]
}
# 这张图只要求得票比例和收入有值,与图 1 的有效样本可不同。
vote_data <- fivethirtyeight::hate_crimes %>%
filter(
!is.na(share_vote_trump),
!is.na(median_house_inc)
) %>%
# 严格超过 50% 才算“过半”;factor() 固定两组的横轴顺序。
mutate(
vote_group = factor(
if_else(share_vote_trump > 0.5, "过半", "未过半"),
levels = c("未过半", "过半")
)
)
# 统计每组记录数,将 n 放在组名下方;标签按组名匹配,不受行顺序影响。
group_n <- vote_data %>% count(vote_group)
group_labels <- setNames(
paste0(group_n$vote_group, "\n(n = ", group_n$n, ")"),
as.character(group_n$vote_group)
)
ggplot(vote_data, aes(x = vote_group, y = median_house_inc, fill = vote_group)) +
# 小提琴显示平滑后的分布形状;trim = FALSE 保留自然收尖的尾部。
# 不设置纵轴 limits,避免把超出观测范围的平滑尾部截掉。
geom_violin(
width = 0.68, trim = FALSE, alpha = 0.66,
colour = "#53656D", linewidth = 0.36
) +
# 箱体显示中位数和四分位距;离群值仍由下一层原始点显示。
geom_boxplot(
width = 0.14, fill = "white", colour = "#304750",
outlier.shape = NA, linewidth = 0.48
) +
# 点仅在水平方向轻微错开,纵轴收入不变;seed 固定重复 Knit 时的位置。
geom_point(
position = position_jitter(width = 0.075, height = 0, seed = 20261002),
shape = 21, fill = "#304750", colour = "white",
stroke = 0.18, size = 1.85, alpha = 0.88
) +
scale_fill_manual(
values = c("未过半" = "#9DC6C7", "过半" = "#E8C3AD"),
guide = "none"
) +
scale_x_discrete(labels = group_labels) +
# 刻度显示为千美元($k);只是改了标签,原始收入数值没有变化。
scale_y_continuous(
labels = scales::label_dollar(scale = 1 / 1000, suffix = "k"),
expand = expansion(mult = c(0.12, 0.13))
) +
labs(
x = "特朗普得票是否过半(2016 年)",
y = "家庭收入中位数(2016 年)"
) +
theme_classic(base_size = 11, base_family = font_family) +
theme(
panel.grid.major.y = element_line(colour = "#E9EDEF", linewidth = 0.3),
axis.line = element_line(colour = "#6F7B80", linewidth = 0.3),
axis.ticks = element_line(colour = "#6F7B80", linewidth = 0.3),
axis.text = element_text(colour = "#33464E"),
axis.title = element_text(colour = "#243941"),
plot.margin = margin(8, 12, 8, 8)
)运行后核对
这一任务不提供完整作图代码。小组在同一份本地 R Markdown
报告中完成三张图:先制作 gapminder
的全年份预期寿命箱线图,再用 R 自带的 airquality
数据制作两张自行选型的图。每张图都应回答一个明确的问题。
图 A 回答:1952—2007 年各大洲国家或地区的预期寿命分布如何变化?
gapminder 的一行是一个国家或地区在某一年的记录。建议先选
lifeExp、continent 和 year
三列:横轴放年份,纵轴放预期寿命(年),按大洲分面。每个“大洲 ×
年份”的箱线图来自该组内各样本国家或地区的观测值;箱线图主要概括中位数、四分位距和箱须外的观测,不能展示全部分布细节。
完成顺序
year 转为因子,或显式设置
group = year。若直接把数值型 year
作为横轴且没有指定分组,geom_boxplot()
默认不会按年份生成多只箱线图。datasets::airquality 收录了 1973 年 5—9
月纽约地区的逐日记录,每行代表一天,包含臭氧、太阳辐射、风速、气温、月份和日期。先用
?airquality
核对变量含义与单位,再从这份数据中提出两个能够用图回答的描述性问题。可从月份间的变化、数值分布或两个变量的关系入手;具体图形由小组决定。
count(continent, year)
检查各组国家数,再决定箱线图是否适合所有分面。factor(year);作图时再检查每个年份是否出现独立箱体。geom_boxplot() 与 facet_wrap()
入手;若补充原始点,应避免让点完全遮住箱体。airquality 可用
summary(datasets::airquality)
检查缺失情况。Month 用数字记录月份;若按日作图,须结合
Month 和
Day,避免把不同月份的同一日混在一起。ggplot 图形末尾加
+ theme(text = element_text(family = "Microsoft YaHei"))(Windows)。macOS
或 Linux 可改用本机已安装的中文字体,例如 PingFang SC 或
Noto Sans CJK SC;单独使用 geom_text()
标注时,也要检查该图层的 family 设置。图 A 参考画法:1952—2007 年各大洲预期寿命的变化。前四洲按年份绘制箱线图;大洋洲每年仅有两个国家,直接显示原始点。
下面展示两种可行的选图方式:图 B 用方格呈现逐日臭氧记录及缺测,图 C 用散点呈现同一天的气温与臭氧记录。小组仍可根据自己提出的问题选择其他图形。
图 B 参考画法:1973 年 5—9 月逐日臭氧浓度。灰色方格表示缺测,右侧标出各月有效天数和当月天数。
图 C 参考画法:当日最高气温与臭氧浓度的散点图。两项指标均有记录的日期共 116 天,颜色区分月份。
以下七题覆盖变量映射、常用图形、时间与坐标尺度、配色和读图判断。每题选一项,再点击“检查答案”查看解析。
© 华东师范大学 社会发展学院 人口研究所 | DAWN 研究组 | yzliu@soci.ecnu.edu.cn
课程负责人:刘贇喆 本章作者:刘贇喆 | 蒋娴静
最后更新:2026年10月14日 构建环境:R version 4.6.1 (2026-06-24 ucrt)