ggplot2:数据操作

ggplot2:数据操作

本文介绍了ggplot2绘图前对数据处理常用的工具。

plyr包简介

ddply(.data, .variable, .fun, …)函数

该函数能够同时在数据的多个子集上作统计汇总。

  • .data:用来作图的数据;
  • .variables:是对数据提取子集的分组变量,形式是.(var1, var2)。为了与图形保持一致,该变量必须包含所有在画图过程中用到的分组变量和分面变量(faceting variable,分面变量用来把数据分割成几个部分,每个部分分别绘制在一个小图里;而分组变量则是将同一个图内的数据分成几部分处理)。
  • .fun是要在各子集上运行的统计汇总函数。这个函数可以返回向量也可以返回数据框。注意这里不要求.fun的返回结果包含分组变量,如果需要的话分组变量会被自动加入到最终的结果里。函数返回结果可以是高度概括的数据集,也许只是一个数,也可以是用修改或拓展之后的数据。

subset()函数

用来对数据取子集的函数,选择数据中的前或后n个或x%的观测值,或是在某个阈值之上或之下的观测值:

1
2
3
4
5
6
7
8
9
## 选取各个颜色里面最小的钻石
library(plyr)
ddply(diamonds, .(color), subset, carat == min(carat))
## 选取最小的两颗钻石
ddply(diamonds, .(color), subset, order(carat) < 2)
## 选取每组里大小为前1%的钻石
ddply(diamonds, .(color), subset, carat > quantile(carat, 0.99))
## 选取所有比平均值大的钻石
ddply(diamonds, .(color), subset, price > mean(price))

transform()函数

用来进行数据变换,与ddply()一起可以计算分组统计量,例如各组的标准差,并且加到原数据上去,例如计算各组的标准差再加到原数据上去。

1
2
3
4
## 把每个颜色组里的钻石的价格标准化
ddply(diamonds, .(color), transform, price = scale(price))
## 减去组均值
ddply(diamonds, .(color), transform, price = price - mean(price))

colwise()函数

colwise()函数用来向量化一个普通函数,也就是说,colwise()能把原本只接受向量输入的函数变成一个可以接受数据框输入的函数。要注意colwise函数返回的是一个新的函数而不是函数运行结果。下面的例子里函数missing()计算向量里缺失值的数目,经过colwise向量化之后,可以应用到数据框上计算数据框各列的缺失值数目。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
> library(ggplot2)
> library(plyr)
> nmissing <- function(x) sum(is.na(x))
> nmissing(msleep$name)
[1] 0
> nmissing(msleep$brainwt)
[1] 27
>
> nmissing_df <- colwise(nmissing)
> nmissing_df(msleep)
name genus vore order conservation sleep_total sleep_rem sleep_cycle
1 0 0 7 0 29 0 22 51
awake brainwt bodywt
1 0 27 0
> colwise(nmissing)(msleep)
name genus vore order conservation sleep_total sleep_rem sleep_cycle
1 0 0 7 0 29 0 22 51
awake brainwt bodywt
1 0 27 0

numcolwise()函数

numclowise()函数是colwise()函数的一个特殊版本,功能类似,但numclowise()只对数值类型的列操作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> msleep2 <- msleep[, -6]
> # 计算中位数
> numcolwise(median)(msleep2, na.rm = T)
sleep_rem sleep_cycle awake brainwt bodywt
1 1.5 0.3333333 13.9 0.0124 1.67
> # 计算分位数(最小值,25%分位数,中位数,75%分位数,最大值)
> numcolwise(quantile)(msleep2, na.rm = T)
sleep_rem sleep_cycle awake brainwt bodywt
1 0.1 0.1166667 4.10 0.00014 0.005
2 0.9 0.1833333 10.25 0.00290 0.174
3 1.5 0.3333333 13.90 0.01240 1.670
4 2.4 0.5791667 16.15 0.12550 41.750
5 6.6 1.5000000 22.10 5.71200 6654.000
> # 计算指定的分位数
> numcolwise(quantile)(msleep2, na.rm = T, probs = c(0.3, 0.4))
sleep_rem sleep_cycle awake brainwt bodywt
1 1.1 0.2000000 11.20 0.0045 0.2984
2 1.4 0.2233333 12.86 0.0066 0.8740

以上的这些函数就可以和ddply一起对数据进行各种分组统计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> # 以vore为分组变量,计算各组各变量的中位数
> ddply(msleep2, .(vore), numcolwise(median), na.rm = T)
vore sleep_rem sleep_cycle awake brainwt bodywt
1 carni 1.95 0.3833333 13.6 0.044500 20.490
2 herbi 0.95 0.2166667 13.7 0.012285 1.225
3 insecti 3.00 0.1666667 5.9 0.001200 0.075
4 omni 1.85 0.5000000 14.1 0.006600 0.950
5 <NA> 2.00 0.1833333 13.4 0.003000 0.122
> # 以vore为分组变量,计算各组各变量的均值
> ddply(msleep2, .(vore), numcolwise(mean), na.rm = T)
vore sleep_rem sleep_cycle awake brainwt bodywt
1 carni 2.290000 0.3733333 13.62632 0.07925556 90.75111
2 herbi 1.366667 0.4180556 14.49062 0.62159750 366.87725
3 insecti 3.525000 0.1611111 9.06000 0.02155000 12.92160
4 omni 1.955556 0.5924242 13.07500 0.14573118 12.71800
5 <NA> 1.880000 0.1833333 13.81429 0.00762600 0.85800

如果以上的函数不够用,读者还可以自己编写函数,只要它能接收、输出数据框即可。例如下面的函数用于计算价格和克拉的秩相关系数以及对数变换后的普通相关系数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
my_summary <- function(df){
with(df, data.frame(
pc_cor = cor(price, carat, method = "spearman"),
lpc_cor = cor(log(price), log(carat))
))
}

# 结果
> ddply(diamonds, .(cut), my_summary)
cut pc_cor lpc_cor
1 Fair 0.9056144 0.9085131
2 Good 0.9599684 0.9687510
3 Very Good 0.9688534 0.9716746
4 Premium 0.9605332 0.9697578
5 Ideal 0.9537275 0.9661884

> ddply(diamonds, .(color), my_summary)
color pc_cor lpc_cor
1 D 0.9561208 0.9606617
2 E 0.9600994 0.9643845
3 F 0.9641572 0.9623876
4 G 0.9633681 0.9696785
5 H 0.9730390 0.9801569
6 I 0.9834392 0.9865118
7 J 0.9846710 0.9879449

拟合多个模型

1
2
3
4
5
6
7
8
9
10
11
12
theme_set(theme_bw(base_size = 30, base_family = "STSong"))
p1 <- qplot(carat, price, data = diamonds, geom = "smooth", colour = color)
dense <- subset(diamonds, carat < 2)
p2 <- qplot(carat, price, data = dense, geom = "smooth", colour = color, fullrange = T)
library(grid)
png("20181021a1.png", width = 1800, height = 900)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
vp <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(p1, vp = vp(1, 1))
print(p2, vp = vp(1, 2))
dev.off()

那么如何使用plyr包生成一个一模一样的图形呢?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(mgcv)
smooth <- function(df){
mod <- gam(price ~ s(carat, bs = "cs"), data = df)
grid <- data.frame(carat = seq(0.2, 2, length = 50))
pred <- predict(mod, grid, se = T)
grid$price <- pred$fit
grid$se <- pred$se.fit
grid
}
(smoothes <- ddply(dense, .(color), smooth))
(q1 <- qplot(carat, price, data = smoothes, colour = color, geom = "line"))
(q2 <- qplot(carat, price, data = smoothes, colour = color, geom = "smooth", ymax = price + 1.96 * se, ymin = price - 1.96 * se))
png("20181021a2.png", width = 1800, height = 900)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
vp <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(q1, vp = vp(1, 1))
print(q2, vp = vp(1, 2))
dev.off()

如果分组变量也被纳入模型的话:

1
2
3
4
5
6
7
mod <- gam(price ~ s(carat, bs = "cs") + color, data = dense)
grid <- with(diamonds, expand.grid(
carat = seq(0.2, 2, length = 50),
color = levels(color)
))
grid$pred = predict(mod, grid)
qplot(carat, pred, data = grid, colour = color, geom = "line")

数据的长宽转换

melt()函数有三个参数:

  1. data:待变形的元数据;
  2. id.vars:依旧放在原列、位置保持不变的变量。id.vars通常是离散的并且预先给定的;
  3. measure.vars:需要被放进同一列的变量。不同的变量放在同一列后,根据原变量名进行分组,这样这些变量就可以同时画在一张图上了。

多重时间序列

1
2
3
4
ggplot(economics, aes(date)) +
geom_line(aes(y = unemploy, colour = "unemploy")) +
geom_line(aes(y = uempmed, colour = "uempmed")) +
scale_colour_hue("variable")


1
2
3
require(reshape2)
emp <- melt(economics, id = "date", measure = c("unemploy", "uempmed"))
qplot(date, value, data = emp, geom = "line", colour = variable)

由于两个变量取值差异太大,所以uempmed变量变成了一条水平线。ggplot2不允许画带两个不同坐标轴的图,因为这样的图有误导性。ggplot2有两个直观的改进方法:把数据调到相同的范围或使用只有标度的分面:

1
2
3
4
5
6
range01 <- function(x){
rng <- range(x, na.rm = T)
(x - rng[1]) / diff(rng)
}
emp2 <- ddply(emp, .(variable), transform, value = range01(value))
qplot(date, value, data = emp2, geom = "line", colour = variable, linetype = variable)


1
qplot(date, value, data= emp, geom = "line") + facet_grid(variable ~ ., scales = "free_y")

平行坐标图

平行坐标图就是以variable变量为x轴表示变量名,以value为y轴表示变量取值。此外,还需要一个分组变量来把各个观测分组(所以每个观测分别对应一条线)。可以直接使用数据库的rownames作为分组变量,但要注意不要让它跟现有的变量名弄混了。我们可以为它取个特殊的名字,比如.row。做好这些准备工作,剩下的就容易了。

1
2
3
4
5
6
7
8
9
10
library(ggplot2movies)
popular <- subset(movies, votes > 1e4)
ratings <- popular[, 7:16]
ratings$.row <- rownames(ratings)
molten <- melt(ratings, id = ".row")

pcp <- ggplot(molten, aes(variable, value, group = .row))
pcp + geom_line(colour = "black", alpha = 1/10)
jitter <- position_jitter(width = 0.25, height = 2.5)
pcp + geom_line(colour = "black", alpha = 1/10, position = "jitter")


为了更清楚地观察电影得分地规律,我们把电影进行聚类,使投票模式相近的被分为一类。下面的代码使用K均值聚类把所有电影分为6类。把平均分最低的记为一类,平均得分最高的记为第六类。

1
2
3
4
5
6
cl <- kmeans(rating[1:10], 6)
ratings$cluster <- reorder(factor(cl$cluster), popular$rating)
levels(ratings$cluster) <- seq_along(levels(ratings$cluster))
molten <- melt(ratings, id = c(".row", "cluster"))
pcp_cl <- ggplot(molten, aes(variable, value, group = .row, colour = cluster))
pcp_cl + geom_line(position = "jitter", alpha = 1/6)


1
pcp_cl + stat_summary(aes(group = cluster), fun.y = mean, geom = "line")

这些图可以很好地看出类之间的差异,但是不能看出聚类的结果好不好。
下图用faceting把每类数据单独画在单独的图里。这张图说明类内的差异还是很大的,说明需要增加类的个数:

1
pcp_cl + geom_line(position = "jitter", colour = "black", alpha = 1/6) + facet_wrap( ~ cluster)

ggplot()方法

ggplot()是个泛型函数,能根据不同的数据调用不同的方法。最常见的数据类型是数据框,所以本章一直以数据框为例。ggplot2()函数更base和lattice系统一样,也可以接收其它数据类型。但是具体的实现方法上有着本质的不同:ggplot2并不像它们一样直接输入数据、输出图形,而只提供作图需要的工具。

ggplot2需要借助fortify()来完成整个过程。fortify()能接收一个模型或对象,以及一个可选的原始数据作为第二参数,然后把它变成能输入到ggplot2的形式,比如数据框。

本节介绍fortify()的原理以及如何用它生成符合ggplot2理念的新方法。ggplot2最重要的理念之一就是数据变形和图形展示尽可能分开进行。这样做可以让函数的使用者不必局限在函数作者设想的某个特殊场合里。比如某些情况下数据变形过程不适用来,图形展示过程还是可以继续使用。

qplot()函数只是对ggplot()的简单封装,所以qplot()也可以接收不同类型的输入。

线性模型

ggplot2目前只有一种适合线性模型的fortify方法。

ggplot2把数据整理和展示完全分离开了。fortify()方法负责数据变形,然后ggplot2负责绘图。fortify()把下表中由模型产生的新变量加到原数据集里。这些变量都是plot.lm()画图中需要产生的。为了避免与数据集中的变量冲突,我们在变量名的前面加上了.

变量 描述
.cooksd Cook距离
.fitted 拟合值
.hat 协方差矩阵估计值的对角线取值
.resid 残差
.sigma 当相应的观测值从模型中被剔除后残差标准差的估计值
.stdresid 标准化残差
1
qplot(displ, cty, data = mpg) + geom_smooth(method = "lm")


1
2
3
4
5
6
mod <- lm(cty ~ displ, data = mpg)
(basic <- ggplot(mod, aes(.fitted, .resid)) +
geom_hline(yintercept = 0, colour = "grey50", size = 0.5) +
geom_point() +
geom_smooth(size = 0.5, se = F))
`


1
(full <- basic %+% fortify(mod, mpg))


1
full + aes(colour = factor(cyl))


1
full + aes(displ, colour = factor(cyl))

编写自己的方法

编写自己的fortify()方法需要知道哪些变量对模型诊断有用以及如何把它们返回给使用者。目前对线性模型的处理方法是直接把它们加到原数据框上,但是有些时候这么做并不好。因为也许用户希望返回一个由数据框组成的列表,每个数据框包含不同的变量,对应不同层次的整合。

fortify()并非一定要和模型返回的对象一起使用。下例就是如何编写一个新的fortify()方法,把一个已有的图片添加到图中。EBImage包负责把图像导入R中,fortify()负责把它变成ggplot2能处理的格式(数据框)。

下面的代码可以绘制我的照片:

1
2
3
4
5
6
7
8
9
10
11
12
13
source("http://bioconductor.org/biocLite.R")
biocLite("EBImage")
library(EBImage)
fortify.Image <- function(model, data, ...){
colours <- channel(model, "x11")
colours <- colours[, rev(seq_len(ncol(colours)))]
melt(colours, c("x", "y"))
}
library(ggplot2)
library(reshape2)
library(ggthemes)
img <- readImage("http://www.czxa.top/MyResume/tex/me.jpg")
qplot(x, y, data = img, fill = value, geom = "tile") + scale_fill_identity() + coord_equal() + theme_map()

# R

评论

程振兴

程振兴 @czxa.top
截止今天,我已经在本博客上写了658.8k个字了!

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×