《ggplot2：数据分析与图形艺术》--学习笔记9

数据操作

plyr包

ddply()函数能够同时在数据的多个子集上作统计汇总。该函数能够根据行的取值，把数据框分解成几个子集，分别把各个子集输入某个函数，最后把结果综合在一个数据框内。基本用法是ddply(.data, .variables, .fun, ...),其中“

• .data 是用来作图的数据。
• .variables是对数据取子集的分组变量，形式是.(var1,var2)。
• .fun是要在各子集上运行的统计汇总函数。
• .subset()用来对数据取子集的函数，选择数据中前（或后）n个（或x%的）观测值，或是在某个阈值之上或下的观测值。

library(ggplot2)
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()是用来进行数据变换的函数。

## 把每个颜色组里钻石的价格标准化，使其均值为0，方差为1
ddply(diamonds,.(color),transform, price = scale(price))
## 减去组均值
ddply(diamonds,.(color),transform,price=price-mean(price))

• colwise()用来向量化一个普通函数，即能把原本只接受向量输入的函数变成可以接受数据框输入的函数。下面例子里函数nmissing()计算向量里缺失值的数目，用colwise()向量化后，可以应用到数据框上，计算数据框各列的缺失值数目。

> nmissing <- function(x) sum(is.na(x))
> nmissing(msleep$name)
[1] 0
> nmissing(msleep$bodywt)
[1] 0
> nmissing(msleep$brainwt)
[1] 27
> nmissing_df <- colwise(nmissing)
> nmissing_df(msleep)
  name genus vore order conservation sleep_total sleep_rem
1    0     0    7     0           29           0        22
  sleep_cycle awake brainwt bodywt
1          51     0      27      0
> ## 更便捷的方法
> colwise(nmissing)(msleep)
  name genus vore order conservation sleep_total sleep_rem
1    0     0    7     0           29           0        22
  sleep_cycle awake brainwt bodywt
1          51     0      27      0

numcolwise()是colwise()的一个特殊版本，功能类似，但只对数值类型的列操作。

> 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
> 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,probs=c(0.25,0.75),na.rm=T)
  sleep_rem sleep_cycle awake brainwt bodywt
1       0.9   0.1833333 10.25  0.0029  0.174
2       2.4   0.5791667 16.15  0.1255 41.750

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

> 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

• 如果以上提供的函数不够用，可以编写自己的函数，只要它也能接收、输出数据框就可以。

> 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

注意，上面例子里的函数my_summary()里没有输出分组变量cut和color，但是最终输出里却自动包含了这些分组变量。

拟合多个模型

怎么用plyr完成ggplot2内置的统计功能？

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)

如何用plyr生成一模一样的图呢？首先赌侠stat_smooth()的文档，看看它内部用的是什么模型：对于大数据用的是gam(y ~ s(x,bs = “cs”))。为了得到与stat_smooth同样的结果，先拟合这个模型，然后在一系列等距的点上求模型的预测值。

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)
qplot(carat,price,data = smoothes,colour = color,geom = "line")
qplot(carat,price,data = smoothes,colour = color,geom = "smooth",
      ymax = price+2*se,ymin=price-2*se)

如果分组变量也被纳入了平滑模型的话，那么用plyr完成这个过程可以更灵活些。

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")

把数据的“宽”变为“长”

ggplot2进行数据分组时必须根据行，而不能根据列。例如可以把钻石根据颜色分组，却不能把钻石的克拉数和价格变量分成两组。
本节将会用到reshape2包的melt()函数和cast()函数。
melt()函数有三个参数：

• data: 待变形的数据。
• id.vars:依旧放在列上、位置保持不变的变量。即数据表的主键。
• measure.vars:需要被放进同一列的变量。不同的变量放在同一列后，根据原变量名来分组。

多重时间序列

时间序列值存在value变量里，然后用vairable变量来区分。

ggplot(economics,aes(date))+
  geom_line(aes(y=unemploy,colour = "unemploy")) +
  geom_line(aes(y=uempmed,colour = "uempemed")) +
  scale_colour_hue("variable")

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

有一个问题：两个变量取值差异太大。ggplot2不允许画带两个不同坐标轴的图，它认为这样的图具有误导性，并给出两个比较直观的改进方法：把数据标度调整到相同的范围（比如：归一化）或使用自由标度的分面。

range01 <- function(x){
  rng <- range(x,na.rm = T)
  (x - rng[1]) / diff(rng)
}
# 将数据框，根据variable分类，对value做归一化变形
emp2 <- ddply(emp, .(variable),transform,value = range01(value))
qplot(date,value,data = emp2,geom = "line",colour = variable,linetype = variable)

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

平行坐标图

即，把原数据化为“长”数据，然后画出平行坐标图。就是以variable变量为x轴表示变量名，以value为y轴表示变量取值。此外，还需要一个分组变量来把各个观测分组。

#取经济数据子集
emp3 <- subset(economics,date>as.Date("1967-08-3"))
#去掉不需要的时间标签
emp3 <- emp3[,-1]
#用ggplot2中的rownames函数为每行加数字标签
emp3$.row <- rownames(emp3)
#对子集中的变量做归一化
emp3$pce <- range01(emp3$pce)
emp3$pop <- range01(emp3$pop)
emp3$psavert <- range01(emp3$psavert)
emp3$unemploy <- range01(emp3$unemploy)
emp3$uempmed <- range01(emp3$uempmed)
#按照rownames为主键融合
molten <- melt(emp3,id=".row")
ggplot(molten,aes(variable,value,group=.row))+geom_line()

透明化

ggplot(molten,aes(variable,value,group=.row))+geom_line(colour="black",alpha=1/20)

增加扰动

ggplot(molten,aes(variable,value,group=.row))+
geom_line(position = 
position_jitter(width = 0.25,height = 2.5))

透明+扰动

ggplot(molten,aes(variable,value,group=.row))+
geom_line(colour="black",alpha=1/20,
position = "jitter")

为了更清楚地观察数据的规律，我们进行聚类，数值相近的分到一类。使用kmeans聚类把所有数值分成6类。

cl <- kmeans(emp3[1:5],6)
#把cl中的cluster按照pce的大小重新整理，并放入emp3中
emp3$cluster <- reorder(factor(cl$cluster),emp3$pce)
#对emp3中的cluster顺序排序
levels(emp3$cluster) <- seq_along(levels(emp3$cluster))
molten <- melt(emp3,id=c(".row","cluster"))
#聚类后所有数据
pcp_cl <- ggplot(molten,aes(variable,value,group=.row,colour = cluster))
pcp_cl + geom_line(position = "jitter",alpha=1/2)

聚类后，每类的均值

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

分面查看每个聚类组内差异不明显，表明聚类效果很好。

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

ggplot()方法

线性模型

fortify方法
fortify()把由模型产生的新变量极大原数据集中。为了避免与数据集中的变量名冲突，我们在变量名的前面都加了“.”

变量	描述
.cooksd	Cook距离
.fitted	拟合值
.hat	帽子矩阵的对角线取值
.resid	残差
.sigma	当相应的观测从模型中被剔除后残差标准差的估计值
.stdresid	标准化残差

一个简单的模型

qplot(displ,cty,data = mpg) + geom_smooth(method = "lm")

用标准化残差代替未标准化残差，让点的大小代表Cook距离。

mod <- lm(cty ~ displ,data = mpg)
p1<-basic <- ggplot(mod,aes(.fitted,.resid))+
  geom_hline(yintercept = 0,colour = "grey50",size=0.5)+
  geom_point()+
  geom_smooth(size=0.5,se=F)
p2<-basic + aes(y = .stdresid)
p3<-basic + aes(size = .cooksd) + scale_size_area("Cook's distance")

（左）基本拟合值-残差图。（中）标准化残差。（右）点大小与Cook距离成比例。

还可以继续增补数据集，把原数据中没有被放入回归模型的变量也添加到图中。这样有助于发现哪些变量有助于改进模型。

full <- basic %+% fortify(mod,mpg)
full + aes(colour = factor(cyl))
full + aes(displ,colour = factor(cyl))

增加了气缸数量变量cyl后，可以看到只要给定了cyl，displ与cty之间其实是呈线性关系的。