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

最新推荐文章于 2023-10-31 16:42:43 发布

Q_M_Y_Y

最新推荐文章于 2023-10-31 16:42:43 发布

阅读量547

点赞数

本文链接：https://blog.csdn.net/zhenxiaojie/article/details/106888776

版权

工具箱

图层叠加的总体策略

图层有3种用途

展示数据本身
展示数据的统计摘要
添加额外的元数据（metadata）、上下文信息和注解

基本图形类型

使用不同的几何对象展示同样数据的效果。

library(ggplot2)
library(grid)
df <- data.frame(
  x = c(3,1,5),
  y = c(2,4,6),
  label = c('a','b','c')
)
center <- theme(plot.title = element_text(hjust = 0.5) )
p <- ggplot(df,aes(x,y)) + xlab(NULL) + ylab(NULL)
p1 <- p + geom_point() + labs(title = 'geom_point') + center
p2 <- p + geom_bar(stat = 'identity') + 
  labs(title = 'geom_bar(stat = \'idensity\')') + center
p3 <- p + geom_line() + labs(title = 'geom_line') + center
p4 <- p + geom_area() + labs(title = 'geom_area') + center
p5 <- p + geom_path() + labs(title = 'geom_path') + center
p6 <- p + geom_text(aes(label = label)) + labs(title = 'geom_text') + center
p7 <- p + geom_tile() + labs(title = 'geom_tile') + center
p8 <- p + geom_polygon() + labs(title = 'geom_ploygon') + center

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,4)))
print(p1,vp = viewport(layout.pos.row = 1,layout.pos.col = 1))
print(p2,vp = viewport(layout.pos.row = 1,layout.pos.col = 2))
print(p3,vp = viewport(layout.pos.row = 1,layout.pos.col = 3))
print(p4,vp = viewport(layout.pos.row = 1,layout.pos.col = 4))
print(p5,vp = viewport(layout.pos.row = 2,layout.pos.col = 1))
print(p6,vp = viewport(layout.pos.row = 2,layout.pos.col = 2))
print(p7,vp = viewport(layout.pos.row = 2,layout.pos.col = 3))
print(p8,vp = viewport(layout.pos.row = 2,layout.pos.col = 4))

在这里插入图片描述

展示数据分布

永远不要指望依靠默认的参数就能对某个具体的分布获得一个表现力强的图形。

depth_dist <- ggplot(diamonds, aes(depth)) + xlim(58,68)
depth_dist + geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  facet_grid(cut ~ .)
depth_dist + geom_histogram(aes(fill = cut), binwidth = 0.1, position = "fill")
depth_dist + geom_freqpoly(aes(y = ..density.., colour = cut), binwidth = 0.1)

直方图
条件密度图
频率多边形图
针对类别型或连续型变量取条件所得的箱线图。

library(plyr)
qplot(cut, depth,data = diamonds, geom = "boxplot")
qplot(carat,depth,data = diamonds, geom = "boxplot",
      group = round_any(carat,0.1,floor),xlim = c(0,3))

在这里插入图片描述
密度图实际上就是直方图的平滑化版本。

qplot(depth,data = diamonds,geom = "density",xlim = c(54,70),
      fill = cut, alpha = I(0.2))

在这里插入图片描述

处理遮盖绘制问题

通过绘制中空，或者更小的点来缓解：
geom_point(shape = 1) ## 中空的点
geom_point(shape = '.') ## 点的大小为像素级
通过透明度调节：
alpha = 1/3
最小的透明度为1/256
通过增加数据扰动：
position_jitter(width = 0.4) ## 增加的默认扰动是数据分辨率的40%
将点分箱并统计每个箱中点的数量
stat_bin2d(binwidth = c(0.02, 200)) ##正方形显示，x轴bin宽度0.02，y轴bin宽度200
stat_binhex(binwidth = c(0.02, 200)) ##六边形显示，同上
使用stat_density2d作二维密度估计

d <- ggplot(diamonds,aes(carat, price)) + xlim(1,3) +
  theme(legend.position = "none")
d + stat_density2d(geom = "point",aes(size = ..density..),contour = F) + scale_size_area()

在这里插入图片描述

曲面图

暂不支持

绘制地图

书中不支持中国地图，CSDN上找了一篇不错的博文可以学习下：“用R语言绘制动态地图，代码奉上!（REmap包详解）”

书中引用的一篇参考文献《用R软件绘制中国分省市地图》

library(maps)
library(mapdata)
par(mar=c(0,0,0,0))
map("china")

在这里插入图片描述

首先，从这里下载中国地图的 GIS 数据，这是一个压缩包，完全解压后包含三个文件（bou2_4p.dbf、bou2_4p.shp 和 bou2_4p.shx），将这三个文件解压到同一个目录下，并在 R 中设好相应的工作空间，然后安装maptools包，运行如下程序：

library(maptools)
x =  readShapePoly('bou2_4p.shp') #下文中会继续用到x这个变量，
                                  #如果你用的是其它的名称，
                                  #请在下文的程序中也进行相应的改动。
par(mar=c(0,0,0,0))
plot(x)

在这里插入图片描述

这时一张完整的中国地图就已经画好了。但是在实际使用的过程中，我们往往会根据自己的需要对地图中的某些省份着以特定的颜色，这时就可以通过调节plot命令中的fg参数来予以实现。然而为了清楚地说明这部分的内容，我需要插播一段 R 绘制地图的原理。

====================== 传说中的分割线 =====================

在绘制地图时，每一个省市自治区或者岛屿都是用一个多边形来表示的。之前的 GIS 数据，其实就是提供了每一个行政区其多边形逐点的坐标，然后 R 软件通过顺次连接这些坐标，就绘制出了一个多边形区域。在上面的数据中，一共包含了 925 个多边形的信息，之所以有这么多是因为一些省份有很多小的附属岛屿。在这 925 个多边形中，每一个都对应一个唯一的 ID，编号分别从 1 到 925。

====================== 传说中的分割线 =====================

回到刚才的话题，plot命令中的fg参数在本例中应该是一个长度为 925 的向量，其第 i 个分量的取值就代表了地图中第 i 个多边形的颜色。一个简单的尝试是运行下面这个命令看看效果：

par(mar=c(0,0,0,0))
plot(x,col=gray(924:0/924));

在这里插入图片描述

在 R 中输入 “x@data”，会得到一个 925 行 7 列的数据框，这其实是 bou24p.dbf 这个文件中存储的信息，之前的readShapePoly()函数虽然读取的是 bou24p.shp 文件，但在默认情况下会把 dbf 文件的信息也放到变量之中。对于这个数据框，其行名就是每一个区域的 ID 编号，第一列和第二列分别是面积和周长，最后一列是该区域所属的行政区名，其它的列应该也是一些编号性质的变量。于是，通过查找相应的行政区对应的行名，就可以对col参数进行赋值了。下面的一个函数，用来生成所需的col向量：

getColor = function(mapdata, provname, provcol, othercol){
	f = function(x, y) ifelse(x %in% y, which(y == x), 0);#如果x的元素在y中存在为真，否则为假，真则返回x的元素在y中的位置，否则位置为0
	colIndex = sapply(mapdata@data$NAME, f, provname);#应用方程f
	col = c(othercol, provcol)[colIndex + 1];#给925个元素按照颜色等级分配颜色
	return(col);
}

其中mapdata是存放地图数据的变量，在上面的例子中就是 x，provname是需要改变颜色的地区的名称，provcol是对应于provname的代表颜色的向量（名称和数字均可），othercol是其它地区的颜色。举例如下：

provname = c("北京市", "天津市", "上海市", "重庆市")
provcol = c("red", "green", "yellow", "purple")
par(mar=c(0,0,0,0))
plot(x, col = getColor(x, provname, provcol, "white"))

在这里插入图片描述

注意provname一定要写地区的全称，写法可以参照下面这条命令生成的向量：

as.character(na.omit(unique(x@data$NAME)))

由此生成的向量有 33 个元素，少了澳门特别行政区，这是这个数据中的一块瑕疵。在x@data的第 899 行有一个NA，不知道它代表的是否就是澳门。

利用类似的方法就可以根据自己的需要对不同的区域进行着色，下面再举一例。从国家统计局获取 2007 年我国各地区的人口数据，然后根据人口的多少对各省份进行着色。程序如下：

provname = c("北京市", "天津市", "河北省", "山西省", "内蒙古自治区",
		"辽宁省", "吉林省", "黑龙江省", "上海市", "江苏省",
		"浙江省", "安徽省", "福建省", "江西省", "山东省",
		"河南省", "湖北省", "湖南省", "广东省",
		"广西壮族自治区", "海南省", "重庆市", "四川省", "贵州省",
		"云南省", "西藏自治区", "陕西省", "甘肃省", "青海省",
		"宁夏回族自治区", "新疆维吾尔自治区", "台湾省",
		"香港特别行政区");
pop = c(1633, 1115, 6943, 3393, 2405, 4298, 2730, 3824, 1858, 7625,
		5060, 6118, 3581, 4368, 9367, 9360, 5699, 6355, 9449,
		4768, 845, 2816, 8127, 3762, 4514, 284, 3748, 2617,
		552, 610, 2095, 2296, 693);
provcol = rgb(red = 1 - pop/max(pop)/2, green = 1-pop/max(pop)/2, blue = 0);
plot(x, col = getColor(x, provname, provcol, "white"), xlab = "", ylab = "");

其中颜色越深的地方代表人口数越多，反之为人口数越少。
在这里插入图片描述

新版本的maptools包的绘图函数已经取消了recs这个参数，现在要实现在地图中只画出这些 ID 所代表的区域，可以在颜色上把不需要的省份变成白色，其中填充色用col参数，边界颜色用border参数。

midchina = c("河南省", "山西省", "湖北省", "安徽省", "湖南省", "江西省")
par(mar=c(0,0,0,0))
plot(x, col = getColor(x, midchina, rep("green", 6),
    "white"), border = "white", xlab = "", ylab = "")

在这里插入图片描述

最后要说的是，在画出的图上仍然可以用points()函数和text()函数加上点和文字，而maptools包中还提供了一个pointLabel()函数，用来解决文本标签的重叠问题。这部分内容请参阅博文：用 R 画中国地图并标注城市位置，以及避免文本标签重叠：maptools 中的 pointLabel()。

在这里插入图片描述

揭示不确定性

用于展示区间的几何对象

变量x类型	仅展示区间	同时展示区间和中间值
连续型	geom_ribbon	geom_smooth(stat=“identity”)
离散型	geom_errorbar geom_linerange	geom_crossbar geom_pointrange

library(ggplot2)
#小于2.5克拉，且随机二项分布（生成钻石长度的0和1的二值分布，每次抽样1的概率0.2）等于1
subset(diamonds,carat<2.5 & rbinom(nrow(diamonds),1,0.2)==1) 
d$lcarat <- log10(d$carat)
d$lprice <- log10(d$price)
#剔除整体的线性趋势
detrend <- lm(lprice ~ lcarat, data = d)
d$lprice2 <- resid(detrend) #提取线性拟合残差

mod <- lm(lprice2 ~ lcarat*color,data = d)#新的线性方程
library(effects)
#以数据框形式展示拟合模型的边际效应结果
effectdf <- function(...){
  suppressWarnings(as.data.frame(effect(...)))
}
#mod模型中color的边际效应
color <- effectdf("color",mod) 
#mod模型中两个参数的边际效应
both1 <- effectdf("lcarat:color",mod)

carat <- effectdf("lcarat",mod,xlevels=50)#将克拉分为50等份
both2 <- effectdf("lcarat:color",mod,xlevels=3)#将每种颜色的克拉分为3个等级

##图5.14：
p1 <- qplot(lcarat,lprice,data = d,colour=color)  #线性拟合
p2 <- qplot(lcarat,lprice2,data=d,colour = color) #拟合残差

library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1,2)))
print(p1,vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
print(p2,vp=viewport(layout.pos.row = 1,layout.pos.col = 2))

## 图5.15：
fplot <- ggplot(mapping = aes(y=fit,ymin=lower,ymax=upper))+
  ylim(range(both2$lower,both2$upper))

p1 <- fplot %+% color + aes(x=color) +geom_point()+geom_errorbar()
p2 <- fplot %+% both2 + aes(x=color, colour = lcarat, group = interaction(color,lcarat)) + 
  geom_errorbar()+geom_line(aes(group=lcarat))+
  scale_colour_gradient()

grid.newpage()
pushViewport(viewport(layout = grid.layout(1,2)))
print(p1,vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
print(p2,vp=viewport(layout.pos.row = 1,layout.pos.col = 2))

## 图5.16：
fplot %+% carat + aes(x=lcarat) + geom_smooth(stat = "identity")

ends <- subset(both1,lcarat==max(lcarat))
fplot %+% both1 + aes(x=lcarat,colour = color)+
  geom_smooth(stat = "identity")+
  scale_colour_hue()+theme(legend.position = "none")+
  geom_text(aes(label=color,x=lcarat+0.02),ends)

在这里插入图片描述 误差棒展示了95%的置信区间

图5.16 模型估计结果中变量carat的不确定性。左carat的边际效应，右针对color的不同水平，变量caret的条件效应。误差带显示了95%的逐点置信区间。

统计摘要

单独的摘要计算函数

midm <- function(x) mean(x,trim = 0.5)

ggplot(diamonds,aes(cut,carat))+
  stat_summary(aes(colour = "trimmed"),fun = midm, geom = "point")+
  stat_summary(aes(colour = "raw"),fun = mean,geom = "point")+
  scale_colour_hue("Mean")

在这里插入图片描述

添加图形注解

unemp <- qplot(date,unemploy,data = economics,geom = "line",xlab = "",ylab = "No. unemployed(1000s)")

在这里插入图片描述

presidential <- presidential[-(1:3), ]

yrng <- range(economics$unemploy)
xrng <- range(economics$date)
unemq + geom_vline(aes(xintercept = as.numeric(start)),data = presidential)

在这里插入图片描述

library(scales)
unemp + geom_rect(aes(NULL,NULL,xmin=start,xmax=end,fill=party),ymin=yrng[1],ymax=yrng[2],
                  data = presidential,alpha=0.2)+scale_fill_manual(values = c("blue","red"))

在这里插入图片描述

last_plot() + geom_text(aes(x=start,y=yrng[1],label=name),
                        data = presidential,size=3,hjust=0,vjust=0)

在这里插入图片描述

caption <- paste(strwrap("Unemployment rates in the US have varied a lot over the years",40),collapse="\n") #字符串按照40分段，中间插入换行后连在一起
unemp + geom_text(aes(x,y,label = caption),data = data.frame(x=xrng[2],y=yrng[2]),
                  hjust = 1,vjust = 1, size = 4) #just=0,0.5,1分别表示左，中，右适应

在这里插入图片描述

highest <- subset(economics,unemploy == max(unemploy))
unemp + geom_point(data = highest, size=3,colour="red",alpha=0.5)

在这里插入图片描述

geom_text：可添加文字描述或为点添加标签。使用取子集的方式，抽取部分观测添加标签可能会非常有用——我们往往希望标注出离群点或其它重要的点。
geom_vline，geom_hline：向图形添加垂直线或水平线。
geom_abline：向图形添加任意斜率和截距的直线。
geom_rect：可强调图形中感兴趣的矩形区域。geom_rect 拥有xmin，xmax，ymin和ymax几种图形属性。
geom_line，geom_path和geom_segment都可以添加直线。所有这些几何对象都有一个arrow参数，可以用来在线上放置一个箭头。我们也可以使用arrow()函数绘制箭头，它拥有angle，length，ends以及type几个参数。

含权数据

权重变量的不同将极大地影响图形内容以及观察结论。用点的大小表示。

p1<-qplot(percwhite,percbelowpoverty,data = midwest)
p2<-qplot(percwhite,percbelowpoverty,data = midwest,
      size = poptotal / 1e6)+
  scale_size_area("Population\n(millions)",breaks = c(0.3,2,4,5))
p3<-qplot(percwhite,percbelowpoverty,data = midwest,size=area) +
  scale_size_area()

grid.newpage()
pushViewport(viewport(layout = grid.layout(1,3)))
print(p1,vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
print(p2,vp=viewport(layout.pos.row = 1,layout.pos.col = 2))
print(p3,vp=viewport(layout.pos.row = 1,layout.pos.col = 3))

在这里插入图片描述 无权重（左图），人口权重（中图），面积权重（右图）

lm_smooth <- geom_smooth(method = lm, size = 1)
p1<-qplot(percwhite,percbelowpoverty,data = midwest)+lm_smooth
p2<-qplot(percwhite,percbelowpoverty,data = midwest,
      weight = popdensity, size = popdensity)+lm_smooth
grid.newpage()
pushViewport(viewport(layout = grid.layout(1,2)))
print(p1,vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
print(p2,vp=viewport(layout.pos.row = 1,layout.pos.col = 2))

在这里插入图片描述 未考虑权重的最优拟合曲线（左图）和以人口数量作为权重的最优拟合曲线（右图）

在我们使用总人口作为权重去修改直方图或密度图的时候，我们的视角将从对郡数量分布的观察转向对人口数量分布的观察。

p1<-qplot(percbelowpoverty,data = midwest,binwidth = 1)
p2<-qplot(percbelowpoverty,data = midwest,weight = poptotal,
      binwidth = 1) + ylab("population")