5.1 统计直方图和核密度估计图
一句话：看单变量的分布：直方图和核密度。
5.1.1 统计直方图
一句话：直方图：数据分组统计频数，最简单常用的分布展示。
统计直方图（histogram），形状类似柱形图却有着与柱形图完全不同的含义。统计直方图涉及统计学的概念，首先要从数据中找出它的最大值和最小值，然后确定一个区间，使其包含全部测量数据，将区间分成若干小区间，统计测量结果出现在各小区间的频数M，以测量数据为横坐标，以频数M为纵坐标，划出各小区间及其对应的频数。在平面直角坐标系中，横轴标出每个组的端点，纵轴表示频数，每个矩形的高代表对应的频数，我们也称这样的统计直方图为频数分布直方图。
所以统计直方图的主要作用如下所示。（1）能够显示各组频数或数量分布的情况。（2）易于显示各组之间频数或数量的差别。
通过统计直方图还可以观察和估计哪些数据比较集中，异常或者孤立的数据分布在何处。统计直方图的基本参数如下所示。（1）组数：在统计数据时，我们把数据按照不同的范围分成几个组，组的个数称为组数。
（2）组距：每一组两个端点的差。（3）频数：分组内的数据元的数量除以组距。
5.1.2 核密度估计图
一句话：核密度估计图：平滑曲线代替直方条，不受分组数影响，分布形状更直观。
核密度估计图（kerneldensityplot）用于显示数据在X轴连续数据段内的分布状况。这种图表是直方图的变种，使用平滑曲线来绘制水平数值，从而得出更平滑的分布。核密度估计图比直方图优胜的地方，在于它们不受所使用分组数量的影响，所以能更好地界定分布形状。
核密度估计（kermeldensityestimation）是在概率论中用来估计未知的密度函数，属于非参数检验方法之一，由Rosenblatt（1955）和EmanuelParzen（1962）[3提出，又名Parzen窗（Parzenwindow）。所谓核密度估计，就是采用平滑的峰值函数（核）来拟合观察到的数据点，从而对真实的概率分布曲线进行模拟。核密度估计，是一种用于估计概率密度函数的非参数方法，xx2x为独立同分布F的n个样本点，设其概率密度函数为f，核密度估计为以下：其中，K0为核函数（非负、积分为1，符合概率密度性质，并且均值为0）。
有很多种核函数，比如高斯函数（gaussian function,f（x)=ae2²，其中a、b和c都为常数），uniform()、triangular()biweight()、triweight()、Epanechnikov()、normal()等。当h>0时，为一个平滑参数，称作带宽（bandwidth）。不同的带宽得到的估计结果差别很大，那么如何选择h？
显然是选择可以使误差最小的。我们用平均积分平方误差（MeanlintergratedSquaredError，MISE）的大小来衡量h的优劣。（b1）单数据系列核密度估计图（b2）多数据系列核密度估计图技能统计直方图和核密度估计图R中的ggplot2包提供了geom_histogram()函数和geom_density()函数，可以分别绘制统计直方图和核密度估计图，
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
df<-read.csv("Hist_Density_Data.csv",stringsAsFactors=FALSE)
#--------------------------------------------(a2) 多数剧系列直方图-----------------------------------
ggplot(df, aes(x=MXSPD, fill=Location))+
geom_histogram(binwidth = 1,alpha=0.55,colour="black",size=0.25)+#, aes(fill = ..count..) )
theme(
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),#,
legend.position=c(0.8,0.8),
legend.background = element_blank()
)
#----------------------------------------(b2)多数剧系列核密度估计图----------------------------------
ggplot(df, aes(x=MXSPD, fill=Location))+
geom_density(alpha=0.55,bw=1,colour="black",size=0.25)+
theme(
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),#,
legend.position=c(0.8,0.8),
legend.background = element_blank()
)
（b2）实现的具体代码如下所示。其中geom_histogram()函数主要由两个参数控制统计分析结果：binwidth（箱形宽度）和bins（箱形总数）：geom_density()函数的主要参数是bw（带宽）和kernel（核函数），核函数默认为高斯核函数“gaussian”，还有其他核函数包括“epanechnikov"，“rectangular"，"triangular"，“biweight”，”cosine”，“optcosine”
峰峦图也可以应用于多数据系列的核密度估计的可视化，如
.png)
.png)
#EasyChartsŶӳƷñؾ
#ʹѧϰϵţEasyCharts
#---------------------------------------ͼ5-1-2ܶȹƷͼ--------------------------------------
library(ggplot2)
library(ggridges)
library(RColorBrewer)
ggplot(lincoln_weather, aes(x = `Mean Temperature [F]`, y = `Month`, fill = ..density..)) +
geom_density_ridges_gradient(scale = 3, rel_min_height = 0.00,size = 0.3) +
scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,'Spectral')))(32))
###------------------------------Զдʵ֣ͼ5-1-2ܶȹƷͼ-------------------------------
library(reshape2)
colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
dt<-lincoln_weather[,c("Month","Mean Temperature [F]")]
splitdata<-split(dt,dt$Month)
xmax<-max(dt$`Mean Temperature [F]`)*1.1
xmin<-min(dt$`Mean Temperature [F]`)*1.1
N<-length(splitdata)
labels_y<-names(splitdata)
mydata<-data.frame(x=numeric(),y=numeric(),variable=numeric()) #յData.Frame
for (i in 1:N){
tempy<-density(splitdata[[i]][2]$`Mean Temperature [F]`,bw = 3.37,from=xmin, to=xmax)
newdata<-data.frame(x=tempy$x,y=tempy$y)
newdata$variable<-i
mydata<-rbind(mydata,newdata)
}
Step<-max(mydata$y)*0.6
mydata$offest<--as.numeric(mydata$variable)*Step
mydata$V1_density_offest<-mydata$y+mydata$offest
p<-ggplot()
for (i in 1:N){
p<-p+ geom_linerange(data=mydata[mydata$variable==i,],aes(x=x,ymin=offest,ymax=V1_density_offest,group=variable,color=y),size =1, alpha =1) +
geom_line(data=mydata[mydata$variable==i,],aes(x=x, y=V1_density_offest),color="black",size=0.5)
}
p+scale_color_gradientn(colours=colormap,name="Density")+
scale_y_continuous(breaks=seq(-Step,-Step*N,-Step),labels=labels_y)+
xlab("Mean Temperature [F]")+
ylab("Month")+
theme_classic()+
theme(
panel.background=element_rect(fill="white",colour=NA),
panel.grid.major.x = element_line(colour = "grey80",size=.25),
panel.grid.major.y = element_line(colour = "grey60",size=.25),
axis.line = element_blank(),
text=element_text(size=15,colour = "black"),
plot.title=element_text(size=15,hjust=.5),
legend.position="right"
)
所示。X轴对应平均温度的数值范围，Y轴对应不同的月份，每个月份的核密度估计数值映射到颜色，这样就可以很好地展示多数据系列的核密度估计结果。0.040.020.00技能核密度估计峰峦图R中的ggridges包提供了geom_density_ridges_gradient()函数，可以结合ggplot2包的ggplot()函数绘制核密度估计峰密图，
的实现代码如下所示。建议将核密度估计峰密图的数值映射到颜色条。1ggridges包的参考手册：https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.htmlgeom_density_ridges_gradient(scale=3,rel_min_height =0.00,size=0.3)+有时候为了更好地发现数据规律或者展示数据分析结果，可以使用二维散点图与统计直方图或核密度估计图的组合图表，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
#-------------------------------------------Method 1: ggpubr包的ggscatterhist()函数------------------------------
library(ggpubr)
N<-300
x1 <- rnorm(mean=1.5, N)
y1 <- rnorm(mean=1.6, N)
x2 <- rnorm(mean=2.5, N)
y2 <- rnorm(mean=2.2, N)
data1 <- data.frame(x=c(x1,x2),y=c(y1,y2))
#(a) 二维散点与统计直方图
ggscatterhist(
data1, x ='x', y = 'y', shape=21,fill="#00AFBB",color = "black",size = 3, alpha = 1,
#palette = c("#00AFBB", "#E7B800", "#FC4E07"),
margin.params = list( fill="#00AFBB",color = "black", size = 0.2,alpha=1),
margin.plot = "histogram",
legend = c(0.8,0.8),
ggtheme = theme_minimal())
N<-200
x1 <- rnorm(mean=1.5, sd=0.5,N)
y1 <- rnorm(mean=2,sd=0.2, N)
x2 <- rnorm(mean=2.5,sd=0.5, N)
y2 <- rnorm(mean=2.5,sd=0.5, N)
x3 <- rnorm(mean=1, sd=0.3,N)
y3 <- rnorm(mean=1.5,sd=0.2, N)
data2 <- data.frame(x=c(x1,x2,x3),y=c(y1,y2,y3),class=rep(c("A","B","C"),each=200))
#(b) 二维散点与核密度估计图
ggscatterhist(
data2, x ='x', y = 'y', #iris
shape=21,color ="black",fill= "class", size =3, alpha = 0.8,
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
margin.plot = "density",
margin.params = list(fill = "class", color = "black", size = 0.2),
legend = c(0.9,0.15),
ggtheme = theme_minimal())
#-----------------------------------Method 2: ggExtra包的ggMarginal()函数------------------------------------
library(ggExtra)
#(a) 二维散点与统计直方图
scatter <- ggplot(data=data1,aes(x=x,y=y)) +
geom_point(shape=21,fill="#00AFBB",color="black",size=3)+
theme_minimal()+
theme(
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=13,face="plain",color="black"),
legend.text= element_text(size=13,face="plain",color="black"),
legend.title=element_text(size=12,face="plain",color="black"),
legend.background=element_blank()
#legend.position = c(0.12,0.88)
)
ggMarginal(scatter,type="histogram",color="black",fill="#00AFBB")
#(b) 二维散点与核密度估计图
scatter <- ggplot(data=data2,aes(x=x,y=y,colour=class,fill=class)) +
geom_point(aes(fill=class),shape=21,size=3)+#,colour="black")+
scale_fill_manual(values= c("#00AFBB", "#E7B800", "#FC4E07"))+
scale_colour_manual(values=c("#00AFBB", "#E7B800", "#FC4E07"))+
theme_minimal()+
theme(
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=13,face="plain",color="black"),
legend.text= element_text(size=13,face="plain",color="black"),
legend.title=element_text(size=12,face="plain",color="black"),
legend.background=element_blank(),
legend.position = c(0.9,0.15)
)
ggMarginal(scatter,type="density",color="black",groupColour = FALSE,groupFill = TRUE)
#-----------------------------------method 3:grid.arrange()函数------------------------------
library(gridExtra)
#(a) 二维散点与统计直方图
# 绘制主图散点图，并将图例去除，这里point层和path层使用了不同的数据集
scatter <- ggplot() +
geom_point(data=data1,aes(x=x,y=y),shape=21,color="black",size=3)+
theme_minimal()
# 绘制上边的直方图，并将各种标注去除
hist_top <- ggplot()+
geom_histogram(aes(data1$x),colour='black',fill='#00AFBB',binwidth = 0.3)+
theme_minimal()+
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank())
# 同样绘制右边的直方图
hist_right <- ggplot()+
geom_histogram(aes(data1$y),colour='black',fill='#00AFBB',binwidth = 0.3)+
theme_minimal()+
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
#axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank())+
coord_flip()
empty <- ggplot() +
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank())
# 要由四个图形组合而成，可以用空白图作为右上角的图形也可以，但为了好玩加上了R的logo，这是一种在ggplot中增加jpeg位图的方法
# logo <- read.jpeg("d:\\Rlogo.jpg")
# empty <- ggplot(data.frame(x=1:10,y=1:10),aes(x,y))+
# annotation_raster(logo,-Inf, Inf, -Inf, Inf)+
# opts(axis.title.x=theme_blank(),
# axis.title.y=theme_blank(),
# axis.text.x=theme_blank(),
# axis.text.y=theme_blank(),
# axis.ticks=theme_blank())
# 最终的组合
grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(4,1), heights=c(1,4))
#(b) 二维散点与核密度估计图
# 绘制主图散点图，并将图例去除，这里point层和path层使用了不同的数据集
scatter <- ggplot() +
geom_point(data=data2,aes(x=x,y=y,fill=class),shape=21,color="black",size=3)+
scale_fill_manual(values= c("#00AFBB", "#E7B800", "#FC4E07"))+
theme_minimal()+
theme(legend.position=c(0.9,0.2))
# 绘制上边的直方图，并将各种标注去除
hist_top <- ggplot()+
geom_density(data=data2,aes(x,fill=class),colour='black',alpha=0.7)+
scale_fill_manual(values= c("#00AFBB", "#E7B800", "#FC4E07"))+
theme_void()+
theme(legend.position="none")
# 同样绘制右边的直方图
hist_right <- ggplot()+
geom_density(data=data2,aes(y,fill=class),colour='black',alpha=0.7)+
scale_fill_manual(values= c("#00AFBB", "#E7B800", "#FC4E07"))+
theme_void()+
coord_flip()+
theme(legend.position="none")
empty <- ggplot() +
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank())
# 要由四个图形组合而成，可以用空白图作为右上角的图形也可以，但为了好玩加上了R的logo，这是一种在ggplot中增加jpeg位图的方法
# logo <- read.jpeg("d:\\Rlogo.jpg")
# empty <- ggplot(data.frame(x=1:10,y=1:10),aes(x,y))+
# annotation_raster(logo,-Inf, Inf, -Inf, Inf)+
# opts(axis.title.x=theme_blank(),
# axis.title.y=theme_blank(),
# axis.text.x=theme_blank(),
# axis.text.y=theme_blank(),
# axis.ticks=theme_blank())
# 最终的组合
grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(4,1), heights=c(1,4))
所示。技能二维散点图与统计直方图组合R中ggpubr包的ggscatterhist()函数（选择"density"参数绘制核密度估计图，选择“histogram"参数绘制统计直方图，选择"boxplot参数绘制箱形图，共三种类型），ggExtra包的ggMarginal()函数（选择"density"参数绘制核密度估计图，选择“histogram”参数绘制统计直方图，选择“boxplot参数绘制箱形图，选择"violin参数绘制小提琴图，共4种类型），gridExtra包的grid.arrange()函数实现ggplot2包绘制的散点图和统计直方图的组合，这三种方法都可以实现二维散点图与统计直方图组合，其中以ggscatterhist()函数最为简单，grid.arrange()函数的可控性最好，也最为复杂。
二维散点与核密度估计图的实现代码如下所示
5.2 数据分布型图表系列
一句话：分布图家族总览：散点/柱形/箱形等多种展示方式。
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20）
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#--------------------------------------图5-2-1四种不同数据分析的分布类图表. (b)核密度估计曲线图.-----------------------------------------------------------
ggplot(mydata, aes(Value,fill=Class))+
geom_density(alpha=1,bw=0.3,colour="black",size=0.25)+
scale_fill_manual(values=brewer.pal(7,"Set2")[c(1,2,4,5)])+
facet_grid(Class~.)+
xlab("X")+
ylab("Desnity")+
theme_light()+
theme(
strip.text = element_text(size=15,color="black"),
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),
legend.position="none"
)
#------------------------------------图5-2-1四种不同数据分析的分布类图表. (a) 统计直方图.-----------------------------------------------
library(reshape2)
type<-as.character(unique(mydata$Class))
step<-0.2
breaks<- seq(min(mydata$Value)-step,max(mydata$Value)+step,step)
mydata1<-data.frame(xvals=numeric(),yvals=numeric(),variable=character()) #创建空的Data.Frame
for (i in 1:length(type)){
x <-mydata[mydata$Class==type[i],2] #rnorm(250 , mean=10 , sd=1)
hg <- hist(x, breaks = breaks , plot = FALSE) # Make histogram data but do not plot
dat <- data.frame(xvals=hg$mids, yvals=hg$counts,variable=rep(type[i],length(hg$mids)))
mydata1 <- rbind(mydata1,dat[dat$yvals>0,])
}
mydata2<-data.frame(xvals=numeric(),value=numeric(),variable=character()) #创建空的Data.Frame
for (i in 1:nrow(mydata1)){
N<-mydata1$yvals[i]
temp<-data.frame(xvals=rep(mydata1$xvals[i],N),value=1:N,variable=rep(mydata1$variable[i],N))
mydata2<-rbind(mydata2,temp)
}
ggplot(mydata2, aes(x=xvals,y=value,fill=variable))+
geom_point(shape=21,size=3,colour="black")+
scale_fill_manual(values=brewer.pal(7,"Set2")[c(1,2,4,5)])+
facet_grid(variable~.)+
xlab("Bins")+
ylab("Count")+
theme_light()+
theme(
strip.text = element_text(size=15,color="black"),
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),
legend.position="none"
)
使用了4种不同数据的分布型数据，每个类别的数据总数分布为100个，其中类别n的数据服从正态分布（normaldistribution：均值μ=3,方差o=1）：类别s的数据为在n数据的基础上右倾斜分布（skew-rightdistribution：Johnson分布的偏斜度2.0和峰度13.1）：类别k的数据在n数据的基础上尖峰态分布（leptikurticdistribution：Johnson分布的偏斜度2.2和峰度20.0）；类别mm为双峰分布（bimodal distribution：两个峰的均值μ1、μ2分别为1.89和3.79，方差o=020.40.20.00.40.00.40.20.40.20.0技能辅助数据的构造使用R自带的rmorm()函数可以构造符合高斯分布的单峰或者多峰数据，使用SuppDists包的rJohnson()函数可以构造符合Johnson分布的数据，然后使用ggplot2包的核密度估计曲线函数geom_density()与分面函数facet_grid()实现如
所示的图表，具体代码如下所示。#生成数据findParams<-function（mu,sigma,skew,kurt）value<-C("JohnsonMomentFitR"as.double(mu),as.double(sigma).as.double(skew).as.double(kurt-3).gamma=double(1).delta =double(1).xi =double(1).lambda =double(1)type=integer(1),PACKAGE=“SuppDists")xi=value$xi,lambda =value$lambda,type=c("SN"“SL",“SU""SB")[value$type])n<-rnorm（100.3.1）#均值为3、标准差为1的正态分布s<-rJohnson（100.findParams（3.1.2.,13.1））#Johnson分布的偏斜度2.0和峰度13.1k<-rJohnson（100,findParams（3,1,2.2,20))#Johnson分布的偏斜度2.2和峰度20.0mm<-rnorm（100,rep（c(2.4）,each=50）*sqrt（0.9），sqrt（0.1））#两个峰的均值μ、μ分别为1.89和3.79，方差a=mydata <-data.frame（Class = factor（rep（c("n","s",“k"."mm”),each =100),.c("n",“s",“K",“mm")).Value =c(n, s,k,#核密度估计曲线图的绘制geom_density(alpha=1,bw=0.3.colour="black"size=0.25)+scale_fill_manual(values=brewer.pal(7,"Set2")[c（1.2,4.5)])+
5.2.1 散点分布图系列
一句话：散点分布图：用散点展示分布，可加误差线/连接线。
散点分布图是指使用散点图的方式展示数据的分布规律，有时可以借助误差线或者连接曲线。
.png)
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20）
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#----------------------------------------------------(a) 散点抖动图---------------------------------------------------------------------------------------
ggplot(mydata, aes(Class, Value))+
geom_jitter(aes(fill = Class),position = position_jitter(0.3),shape=21, size = 2)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
# #------------------------------------------------------(b) 蜂群图---------------------------------------------------------------------------------
library(ggbeeswarm) #library(beeswarm) biocLite(c("beeswarm","ggplot2"))
ggplot(mydata, aes(Class, Value))+
geom_beeswarm(aes(fill = Class),shape=21,colour="black",size=2,cex=2)+
scale_fill_manual(values= c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
xlab("Class")+
ylab("Value")+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#----------------------------------------------------------(c)点阵图----------------------------------------------------------------------------------------
ggplot(mydata, aes(Class, Value))+
geom_dotplot(aes(fill = Class),binaxis='y', stackdir='center', dotsize = 0.6)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------------------------------(e) 带误差线散点与点阵组合图--------------------------------------------
ggplot(mydata, aes(Class, Value,fill = Class))+
geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.6)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
geom_pointrange(stat="summary", fun.data="mean_sdl",fun.args = list(mult=1),
color = "black",size = 1.2)+
geom_point(stat="summary", fun.y="mean",fun.args = list(mult=1),
color = "white",size = 4)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
ggplot(mydata, aes(Class, Value,fill = Class))+
geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.6)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="pointrange", color = "black",size = 1.2)+
stat_summary(fun.y="mean", fun.args = list(mult=1),
geom="point", color = "white",size = 4)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------------------------------(d) 带误差线的散点与抖动图--------------------------------------------
ggplot(mydata, aes(Class, Value))+
geom_jitter(aes(fill = Class),position = position_jitter(0.3),shape=21, size = 2,color="black")+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="pointrange", color = "black",size = 1.2)+
stat_summary(fun.y="mean", fun.args = list(mult=1),
geom="point", color = "white",size = 4)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------------------------------(f)带连接线的带误差线散点图--------------------------------------------
library(ggalt)
library(dplyr)
mydata2 <- mydata %>%
group_by(Class) %>%
summarise(sd = sd(Value),len = mean(Value))
ggplot(mydata2, aes(x = c(1:4), y = len, ymin = len-sd, ymax = len+sd))+
geom_xspline(spline_shape = -0.5,size=1) +
geom_errorbar(colour="black", width=0.2,size=1)+
geom_point(aes(fill = Class),shape=21,size=5,stroke=1)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
scale_y_continuous(breaks=seq(0,7,2),lim=c(0,7.5))+
xlab("Time")+
ylab("Value")+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
所示为6种不同形式的散点分布图。
为抖动散点图（jitterchart），每个类别数据点的Y轴数值保持不变，数据点X轴数值沿着X轴类别标签中心线在一定范围内随机生成，然后绘制成散点图。所以，抖动散点图的主要绘制参数就是数据点的抖动范围。由于随机生成数据点的X轴数值，所以很容易存在数据点重合叠加的情况，不利于观察数据的分布规律。
ggplot2包的geom_jiter()函数可以绘制抖动散点图，其关键参数是position=position_jitter（width=NULL），width表示水平方向左右抖动的范围。
为蜂巢图（beeswarmchart），每个类别数据点沿着X轴类别标签中心线向两侧，同时逐步向上均匀而对称展开，整体较为美观，也方便读者观察数据的分布规律。可以借助ggplot2的拓展包ggbeeswarm中的geombeeswarm()函数，主要参数包括散点的形状（shape）、大小（size）和间隙（cex）。
为点阵图（dotplot），每个类别数据点沿着X轴类别标签中心线向两侧均匀而对称地展开，整体较为美观，很方便读者观察数据的分布规律。ggplot2包的geomdotplot()函数可以绘制点阵图，主要参数包括binwidth（箱形宽度）binaxis（箱形的排布方向）（沿X或Y轴）stackdir（散点的排布方式）（默认为"up"，还有“down"、“center”）、dotsize（散点大小）等。
为抖动散点图+带误差线的散点图，先根据每个类别数据直接绘制散点图，然后添加每个类别数据的均值与误差线（标准差）：average+standarddeviation。如果只使用带误差线的散点图，就无法观察数据的分布情况，所以使用抖动散点图作为背景，可以很好地显示数据分布情况。数据均值与误差线的添加可以使用statsummary()函数实现。
具体地说，即stat_summary（fun.data="mean_sdl"，geom="pointrange"）函数可以绘制带均值点的误差线图。
为点阵图+带误差线的散点图，先根据每个类别的数据直接绘制散点图，然后添加每个类别数据的均值与误差线（标准差）：average+standarddeviation。如果只使用带误差线的散点图，就无法观察数据的分布情况，所以使用点阵图作为背景，可以很好地显示数据分布情况，与5-2-1图（d）表达的信息类似。
（f为带连接线的带误差线散点图，使用曲线连接散点，但是这时的X轴变量为连续型的时间变量，而不是
的类别变量。用曲线连接数据点可以表示数据的变化关系与趋势，与节中的散点曲线图系列基本类似，但此处是添加误差线表示数据的分布情况。我们可以先借助dplyr包的group_by()函数和summarise()函数分组计算不同类别的均值与标准差；然后使用ggplot2包的geom_point()函数和geom_errorbar()函数分别绘制均值点和对应的误差线：最后使用ggalt包的geom_xspline()函数用光滑的曲线连接各点。
技能散点分布图系列
类似，都是带误差线的散点图与分布类散点图的组合，就是使用geom_jitter()函数或者geom_dotplot()函数绘制点阵图或抖动散点图，再添加误差线和均值点。其中图5-2-2（d）的实现代码如下所示。#添加抖动散点geom_jitter(aes（fill =Class).,position=position_jitter(0.3),shape=21, size=2,color=black")+scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+#添加误差线stat_summary(fun.data="mean_sdl",fun.args=list（mult=1).geom=pointrange”,color=“black”size=1.2)+#添加均值散点
5.2.2 柱形分布图系列
一句话：柱形分布图：用柱形展示分布，可叠加散点/误差线。
柱形分布图系列是指使用柱形图的方式展示数据的分布规律，有时可以借助误差线或者散点图。如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20）
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#--------------------------------------------------图5-2-3 柱形分布图系列。(a) 带误差线的柱形图------------------------------------------
ggplot(mydata, aes(Class, Value))+
stat_summary(mapping=aes(fill = Class),fun.y=mean, fun.args = list(mult=1),geom='bar',colour="black",width=.7) +
stat_summary(fun.data = mean_sdl, fun.args = list(mult=1),geom='errorbar', color='black',width=.2) +
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
ylim(0,7.5)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#---------------------------------------图5-2-3 柱形分布图系列(b) 带误差线柱形与抖动图----------------------------------------
ggplot(mydata, aes(Class, Value))+
stat_summary(fun.y=mean, fun.args = list(mult=1),geom='bar',colour="black",fill="white",width=.7) +
stat_summary(fun.data = mean_sdl,fun.args = list(mult=1), geom='errorbar', color='black',width=.2) +
geom_jitter(aes(fill = Class),position = position_jitter(0.2),shape=21, size = 2,alpha=0.9)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
所示。带误差线的柱形图就是使用每个类别的均值作为柱形的高度，再根据每个类别的标准差绘制误差线，如
所示。但是如果只使用
展示数据，那么就会与带误差线的散点图存在同样的问题：无法显示数据的分布情况。
类别nn为双峰分布，但是其与其他三个类别的均值与标准差基本相同，没有较大区别。所以可以在带误差线的柱形图的基础上，添加抖动散点图，这样可以方便观察数据分布规律。技能柱形分布图系列
带误差线柱形与抖动组合图就是在带误差线柱形图的基础上，再使用geom_jitter()函数添加抖动散点图。其中，带误差线柱形图使用stat_summary（fun.y-mean，geom=bar）实现柱形图，而stat_summary（fun.data=mean_sdl,geom=errorbar）实现误差线的绘制。#添加柱形图stat_summary(fun.y=mean,geom='bar',fun.args=list（mult=1),colour="black"fill="white”,width=.7)+#添加误差线stat_summary（fun.data=mean_sdl,fun.args=list（mult=1)geom='errorbar',color=black,width=.2)+#添加抖动散点图geom_jitter(aes(fill =Class).position =position_jitter(0.2),shape=21,size =2.alpha=0.9)+scale_fil_manual(values=c(brewer.pal(7,"Set2")[c(1.2,4,5))+
5.2.3 箱形图系列
一句话：箱形图：最大值/最小值/中位数/四分位——组间对比最常用。
箱形图（boxplot）也称箱须图（box-whiskerplot）、箱线图、盒图，能显示出一组数据的最大值、最小值、中位数，以及上下四分位数，可以用来反映一组或多组连续型定量数据分布的中心位置和散布范围，因形状如箱子而得名。1977年，箱形图首先出现在美国著名数学家JohnW.Tukey的著作ExploratoryData Analysis中3l。它能方便显示数字数据组的四分位数。
从盒子两端延伸出来的线条称为“晶须”（whisker），用来表示上、下四分位数以外的变量。异常值（outlier）有时会以与晶须处于同一水平的单一数据点表示。这种箱形图以垂直或水平的形式出现，如
所示。异常值晶须晶须下极限下四分位数中位数上四分位数上极限其中，四分位数（guartile）是指在统计学中把所有数值由小到大排列并分成四等份，处于三个分割点位置的数值。分位数是将总体的全部数据按大小顺序排列后，处于各等分位置的变量值。
如果将全部数据分成相等的两部分，它就是中位数：如果分成四等分，就是四分位数：八等分就是八分位数等。四分位数也被称为四分位点，它是将全部数据分成相等的四部分，其中每部分包括25%的数据，处在各分位点的数值就是四分位数。四分位数有三个，第一个四分位数就是通常所说的四分位数，也被称为下四分位数，第二个四分位数就是中位数，第三个四分位数称为上四分位数，分别用Q1、Q2、Q3表示。
第一个四分位数（Q1），又被称“较小四分位数”，等于该样本中所有数值由小到大排列后第25%的数字。第二个四分位数（Q2），又被称“中位数”，等于该样本中所有数值由小到大排列后第50%的数字。第三个四分位数（Q3），又被称“较大四分位数”，等于该样本中所有数值由小到大排列后第75%的数字。
第三个四分位数与第一个四分位数的差距又被称为四分位距（InterQuartileRange，IQR），是上四分位值Q3与下四分位值Q1之间的差，即IQR=Q3-Q1。IQR乘以因子0.7413得到标准化四分位距（NormIQR），它是稳健统计技术处理中用于表示数据分散程度的一个量，其值相当于正态分布中的标准偏差（SD）。
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
ggplot(mydata, aes(Class, Value))+
geom_boxplot(aes(fill = Class),notch = FALSE)+
geom_jitter(binaxis = "y", position = position_jitter(0.3),stackdir = "center",dotsize = 0.4)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
所示为箱形图系列。从箱形图可以得出的观察结果主要体现在5个方面：①关键数值，例如平均值、中位数和上下四分位数等：2任何异常值（以及它们的数值）：③数据分布是否对称；4数据分组有多紧密：5数据分布是否出现偏斜（如果是，那么往什么方向偏斜）。箱形图通常用于描述性统计，是以图形方式快速查看一个或多个数据集的好方法。
虽然与直方图或密度图相比似乎有点原始，但它们占用较少空间，当要比较很多组或数据集之间的分布时便相当有用。箱形图在数据显示方面会受到限制，简单的设计往往隐藏了有关数据分布的重要细节，例如在使用箱形图时，我们不能了解数据分布是双模还是多模的。箱形图作为描述统计的工具之一，其功能有独特之处，主要有以下几点。
（1）直观明了地识别批量数据中的异常值。一批数据中的异常值值得关注，忽视异常值的存在是十分危险的，不加剔除地把异常值包括在数据的计算分析过程中，对结果会带来不良影响；重视异常值的出现，分析其产生的原因，常常成为发现问题进而改进决策的契机。箱形图为我们提供了识别异常值的一个标准：异常值被定义为小于Q1-1.5IQR或大于Q3+1.5IQR的值。
虽然这种标准有点任意性，但它来源于经验判断，经验表明它在处理需要特别注意的数据方面表现不错。这与识别异常值的经典方法有些不同。众所周知，基于正态分布的3o法则或z分数方法是以假定数据服从正态分布为前提的，但实际数据往往并不严格服从正态分布。
它们判断异常值的标准是以计算批量数据的均值和标准差为基础的，而均值和标准差的耐抗性极小，异常值本身会对它们产生较大影响，这样产生的异常值个数不会多于总数的0.7%。显然，应用这种方法于非正态分布数据中判断异常值，其有效性是有限的。而箱形图有两方面优势：一方面，其绘制依靠实际数据，不需要事先假定数据服从特定的分布形式，没有对数据做任何限制性要求，它只是真实直观地表现数据形状的本来面貌；另一方面，箱形图判断异常值的标准以四分位数和四分位距为基础，四分位数具有一定的耐抗性，多达25%的数据可以变得任意远而不会很大地干扰四分位数，所以异常值不能对这个标准施加影响，箱形图识别异常值的结果比较客观。
由此可见，箱形图在识别异常值方面有一定的优越性。（2）利用箱形图判断批量数据的偏态和尾重。比较标准正态分布、不同自由度的分布和非对称分布数据的箱形图的特征，可以发现：对于标准正态分布的大样本，只有0.7%的值是异常值，中位数位于上、下四分位数的中央，箱形图的方盒关于中位线对称。
选取不同自由度的t分布的大样本，代表对称重尾分布，当t分布的自由度越小时，尾部越重，就有越大的概率观察到异常值。以卡方分布作为非对称分布的例子进行分析，我们发现当卡方分布的自由度越小时，异常值出现于一侧的概率越大，中位数也越偏离上、下四分位数的中心位置，分布偏态性越强。若异常值集中在较小值一侧，则分布呈现左偏态：若异常值集中在较大值一侧，则分布呈现右偏态。
箱形图可以很好地用于观察数据的分布，但是无法适用于双峰及多峰分布的数据，如
所示类别mm（数据服从双峰分布），可以准确获得数据的分布情况，所以在箱形图的基础上添加抖动散点图或者点阵图，可以方便读者观察原始数据的分布情况，如
所示。技能箱形图系列R中ggplot2包的geom_boxplot()函数可以绘制箱形图，geom_jitter()函数可以绘制抖动散点图，具体代码如下所示。geom_boxplot(aes（fill=Class),notch=FALSE)+geom_jitter（binaxis=“y".position=position_jitter（0.3),stackdir=“center"dotsize=0.4)+scale_fill_manual(values=c（brewer.pal(7,"Set2")[c（1,2,4,5)])+最常用的两种箱形图：可变宽度( variable-width )和带凹槽(notched )的箱形图[32,3],如图 5-2-6(a)和
所示。箱形图的另外一个变量：箱形图的宽度，就是为了解决箱形图每个类别的数据量大小不同的问题[32.33]。
（a)的类别a、b、c和d都服从正态分布，其数据量大小分别为10、100、1000和10000，箱形的宽度依次增加。在
所示带凹槽的箱形图中，中位数的置信区间（confidenceinterval）可以由凹槽对应表示。因此，不考虑数据的分布情况，如果凹槽不重合，表示中位数在95%的置信区间内就可以认为显著不同。技能箱形图系列
可变宽度的带凹槽箱形图可以将geomboxplot()函数的参数notch设置为是否带凹槽（TRUE/FALSE），参数varwidth可以设置为是否根据将箱形宽度映射到箱形宽度（TRUE/FALSE），具体代码如下所示
）能有效地展示数据的分布情况与异常值。但是对于中等数据集（n<1000），对四分位数之外数据的估计可能不可靠，所以箱形图所提供的信息在四分位数之外的情况下是相当模糊的，而对于一个数据集大小为n的高斯样本来说，异常值（outlier）和远外值（far-outvalue）通常小于1034]。而我们希望使用大数据集（n=10，000-100，000）可以提供更加精准的四分位数之外的数据估计，同时可以展示大量的异常值（约0.4+0.007n）。
letter-value箱形图就能满足我们的需求，它不仅能展示四分位数之外的数据分布信息，还能显示异常值的分布情况。letter-value箱形图在箱形图的中值（median（M））和四分位数（fourths（F））的基础上，往两端延伸，增加箱形的个数：1/8eigths（E），1/16sixteenths（D），直到估计误差增大到一定的值。如
所示，一系列的小箱子堆积而成，展示数据的分布情况。但是它与传统箱形图存在一个同样的问题，即无法识别多峰分布的情况[35,36]在
（a)中，类别a、b、c和d都服从正态分布，其数据量大小分别为100、1000、10000和100000。在
中，类别n、s、k和nm服从不同的数据分布，其数据量大小分别为100、1000、10000和100000，其中nm数据服从双峰分布，但是仅仅从图中无法识别，这就是箱形图的局限性所在。对于实验数据的分析与展示，很多人会使用常见的带误差线的柱形图，因为使用Excel就可以直接绘制。但是这样展示数据，信息量是非常低的。
而使用箱形图能够提供更多的数据分布信息，能更好地展示数据（Excel2016版本也提供了箱形图的绘制功能）。在期刊NatureMethods2013年的文章中有100个带误差线的柱形图，而只有20个箱形图，从这里就可以看出来，用箱形图的人远远没有使用带误差线的柱形图的人多。于是自然出版集团（NaturePublishingGroup）写了两篇专栏文章PointsofView:Barcharts and boxplots37Points ofSignificance:Visualizingsamples with boxplots38，并且还发表了一篇文章BoxPlotR：awebtoolforgenerationofboxplots39]，专门对比箱形图与带误差线的柱形图在数据分布展示方面的差异，最后得出的结论是：箱形图能够比带误差线的柱形图更好地展示数据的分布情况。
5.2.4 其他图表
一句话：瓶状图等变体：用核密度曲线替代箱体，看分布形状。
瓶状图（vaseplot）就是使用核密度估计箱形部分的数据，从而得到核密度估计曲线，替代原有的箱形部分，主要用来显示数据的分布形状[4（见
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
#Reference：https://github.com/hadley/boxplots-paper
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#-----------------------------------------------(a) 瓶状图---------------------------------------------------------------------------
source("lvplot/boxplots-vase.r")
#pdf("images/four-vase.pdf", width = 4, height = 4)
par(mar = c(2.1, 2.1, .1, .1))
vase(split(mydata$Value, mydata$Class), bw = rep(0.1, 4))
axis(side = 2)
xlab("Class")
#-------------------------------------------(b)小提琴图---------------------------------------------------------------------------
ggplot(mydata, aes(Class, Value))+
geom_violin(aes(fill = Class),trim = FALSE)+
geom_boxplot(width = 0.2)+
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#-------------------------------------------(c)豆状图---------------------------------------------------------------------------
library(beanplot)
par(mar = c(2.1, 2.1, .1, .1))
beanplot(Value ~Class, data = mydata,col=c("white","black"),xlab ="Class",ylab ="value")
#--------------------------------------------(d) 海盗图-------------------------------------------------------------------------
library(yarrr)
pirateplot(formula = mydata$Value~mydata$Class, data =mydata,
theme.o = 2,
xlab = "Class", ylab = "Value", main = "",
# Choose your color palette, or give common color vector
#pal = color,#"google",
#gl.col = gray(.8),
# Set transparency of the elements:
#bean.b.col="black",
bean.f.col=brewer.pal(7,"Set2")[c(1,2,4,5)],
bar.b.col="black",
#line.o = 0.1,
bar.o = .1,
bean.o = .1,
point.o = .9,
# Shape of point
#point.pch = 2,
#Background color
#back.col = "white",
ylim=c(0,7.5),
gl.col = "white", # gridlines
gl.lwd = c(.5, 0)
) # turn off minor grid lines)
）。绘图时需要设定核密度估计的带宽小提琴图（violinplot）用于显示数据分布及其概率密度（见
）。这种图表结合了箱形图和密度图的特征，主要用来显示数据的分布形状。中间的黑色粗条表示四分位数范围，从其中延伸出的幼细黑线代表95%置信区间，而黑色横线则为中位数4]。
虽然小提琴图可以比箱形图显示更多详情，但它们也可能包含较多干扰信息，而且绘图时需要设定核密度估计的带宽。小提琴图可以使用ggplot2包的geom_violin()函数，主要参数与核密度估计曲线一样，也是bw（带宽）。豆状图（beanplot）是在小提琴密度部分的基础上，用短横条表示每个数据的数值，用长横线表示类别数据的均值（见
（c)）。它看起来就是豌豆，而里面的短横条看起来像里面的种子[42]。豆状图可以使用beanplot包的beanplot()函数实现。
海盗图（pirateplot），中文名是笔者给起的，因为它算是综合了抖动散点图（原始数据）柱形图（均值）小提琴图（核密度估计）和方块—95%的高密度区间（HighDensityInterval，HDI）或置信区间（ConfidenceInterval，CI），虽然能完整地表达数据的所有信息，但是过于复杂43]（见图5-2-8（d））。建议将该图表用于前期数据分布的探索，再具体确定选择合适的图表类型展示数据。海盗图可以使用yarr包的pirateplot()函数实现。
梯度图（gradientplot，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
#Reference：https://github.com/hadley/boxplots-paper
# #--------------------------------------------------数据的准备与导入------------------------------------------------------------------
library(RColorBrewer)
library(ggplot2)
mydata<-read.csv("Norm.csv",stringsAsFactors=FALSE,header=TRUE) #以3为均值，1为标准差的正态分布
mydata$Class<-rep("Class",nrow(mydata))
colnames(mydata)<-c("Value","Class")
color<-brewer.pal(7,"Set2")[c(1,2,4,5)]
#----------------------------------------------------(a)散点抖动图-------------------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_jitter(fill =color[4],position = position_jitter(0.2),shape=21, size = 3)+
scale_y_continuous(breaks=seq(0,6,1))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
axis.title.x = element_blank(),
legend.position="none"
)
p
#-----------------------------------------------------(b) 蜂巢图-----------------------------------------------------------------
library(beeswarm)
class<-mydata$Class
value<-mydata$Value
beeswarm<-beeswarm(value~class, data = mydata,method = 'swarm')[, c(1, 2, 6)]
colnames(beeswarm) <- c("x", "y","Class")
ggplot(beeswarm, aes(x,y)) +
geom_point(fill = color[4],shape=21,colour="black",size=3.5)+
scale_fill_manual(values= c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
scale_y_continuous(breaks=seq(0,6,1))+
xlab("Class")+
ylab("Value")+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#----------------------------------------------(c)点状图--------------------------------------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_dotplot(fill =color[4],binaxis='y', stackdir='center', dotsize = 0.8)+
scale_y_continuous(breaks=seq(0,6,1))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#----------------------------------------------(d)统计直方图--------------------------------------------------------------------
ggplot(mydata, aes(Value, fill=Value))+
geom_histogram(alpha=1,fill=color[4],colour="black",size=0.25)+#binwidth = 1,
coord_flip()+
theme_classic()+
scale_x_continuous(breaks=seq(0,6,1))+
theme(
panel.background = element_rect(color="black"),
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),
legend.position=c(0.8,0.8),
legend.background = element_blank()
)
#----------------------------------------------(e)核密度估计图--------------------------------------------------------------------
ggplot(mydata, aes(Value, fill=Value))+
geom_density(alpha=1,colour="black",size=0.25,fill=color[4])+
coord_flip()+
theme_classic()+
scale_x_continuous(breaks=seq(0,6,1))+
theme(
panel.background = element_rect(color="black"),
text=element_text(size=15,color="black"),
plot.title=element_text(size=15,family="myfont",face="bold.italic",hjust=.5,color="black"),
legend.position=c(0.8,0.8),
legend.background = element_blank()
)
#-------------------------------------------(f)带误差线的散点图-------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_dotplot(fill="white",binaxis='y', stackdir='center', dotsize = 0.8)+
stat_summary(fill = color[4],fun.data="mean_sdl", fun.args = list(mult=1),
geom="pointrange", color = "black",size =2 ,shape=21)+
scale_y_continuous(breaks=seq(0,6,1))+
ylab("Value")+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#-----------------------------------------------(g)带误差线的柱形图----------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
stat_summary(fill =color[4],fun.y=mean, geom='bar',colour="black",width=.7,size=0.5) +
stat_summary(fun.data = mean_sdl, geom='errorbar', color='black',width=.2,size=0.5) +
geom_jitter(fill ="white",position = position_jitter(0.2),shape=21, size = 2,alpha=0.9)+
scale_y_continuous(breaks=seq(0,6,1))+
ylab("Value")+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#------------------------------------------------(h) 梯度图--------------------------------------------------
library(denstrip)
#pdf("images/four-denstrip2.pdf", width = 4, height = 4)
par(mar = c(2.1, 2.1, .1, .1))
plot(c(0.5, 1.5), range(mydata$Value), type = "n", axes = F, xlab = "Class", ylab = "Value")
rect(-10, -10, 10, 10, col = "white")
denstrip(mydata$Value, at = 1,mticks=mean(mydata$Value), hor = F, width = 0.75, bw = 0.2, colmax=brewer.pal(7,"Set2")[c(5)],
colmin="white",mlen=1.1,mcol="black")
box()
axis(2)
axis(1, at = 1, labels = "Class")
#------------------------------------------------(i) 箱型图--------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_boxplot(fill =color[4],notch = FALSE) +
theme_classic()+
scale_y_continuous(breaks=seq(0,6,1))+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#------------------------------------------------(j) 带凹槽的箱型图--------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_boxplot(fill =color[4],notch = TRUE) +
theme_classic()+
scale_y_continuous(breaks=seq(0,6,1))+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#---------------------------------------------------------(k)瓶状图----------------------------------------------------
source("boxplots-vase.r")
par(mar = c(2.1, 2.1, .1, .1))
vase(split(mydata$Value, mydata$Class), bw = 0.15)
axis(side = 2)
xlab("Class")
#------------------------------------------------------ (l)豆状图----------------------------------------------------
library(beanplot)
par(mar = c(2.1, 2.1, 2.1, 2.1))
beanplot(Value ~Class, data = mydata,col=color[4],xlab ="Class",ylab ="value")#,ylim =c(-3,3))
#------------------------------------------------------(m)小提琴图------------------------------------------------
p <- ggplot(mydata, aes(Class, Value))+
geom_violin(fill =color[4],trim = FALSE)+
geom_boxplot(width = 0.2)+
scale_y_continuous(breaks=seq(0,6,1))+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
p
#----------------------------------------------------(n)海盗图--------------------------------------------------------------------
library(yarrr)
head(pirates)
pirateplot(formula = mydata$Value~mydata$Class, data =mydata,
theme.o = 2,
xlab = "", ylab = "Value", main = "",
# Choose your color palette, or give common color vector
pal = color[1],#"google",
#gl.col = color[4],
# Set transparency of the elements:
#bean.b.col="black",
bean.f.col=color[4],
bar.b.col="black",
line.o = 0.9,
bar.o = .4,
bean.o = .1,
point.o = .9,
# Shape of point
#point.pch = 2,
#Background color
#back.col = "white",
ylim=c(0,6),
gl.col = "white", # gridlines
gl.lwd = c(.5, 0)
)
所示），也可以表示数据分布情况，彩色的条带对应数据的核密度估计，黑色长条代表数据的均值或者中位数，表达的数据信息与小提琴图类似44。梯度图可以使用denstrip包的denstrip()函数实现。技能小提琴图
所示的小提琴图也是很常见的图表，可以使用ggplot2包的geom_violin()函数实现。一般我们还可以在小提琴图里添加箱形图，这样能更加全面地展示数据，其核心代码如下所示。geom_violin(aes(fill =Class),trim=FALSE)+scale_fill_manual(values=c(brewer.pal(7.Set2")[c(1,2,4,5)])+云雨图，除以上图表外，推荐大家使用云雨图，云雨图可以看成核密度估计曲线图、箱形图和抖动散点图的组合图表，清晰、完整、美观地展示了所有数据信息，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(grid)
library(RColorBrewer)
library(dplyr)
library(SuppDists) #提供rJohnson()函数
# somewhat hackish solution to:
# https://twitter.com/EamonCaddigan/status/646759751242620928
# based mostly on copy/pasting from ggplot2 geom_violin source:
# https://github.com/hadley/ggplot2/blob/master/R/geom-violin.r
"%||%" <- function(a, b) {
if (!is.null(a)) a else b
}
color<-brewer.pal(7,"Set2")[c(1,2,4,5)]
geom_flat_violin <- function(mapping = NULL, data = NULL, stat = "ydensity",
position = "dodge", trim = TRUE, scale = "area",
show.legend = NA, inherit.aes = TRUE, ...) {
layer(
data = data,
mapping = mapping,
stat = stat,
geom = GeomFlatViolin,
position = position,
show.legend = show.legend,
inherit.aes = inherit.aes,
params = list(
trim = trim,
scale = scale,
...
)
)
}
GeomFlatViolin <-
ggproto("GeomFlatViolin", Geom,
setup_data = function(data, params) {
data$width <- data$width %||%
params$width %||% (resolution(data$x, FALSE) * 0.9)
# ymin, ymax, xmin, and xmax define the bounding rectangle for each group
data %>%
group_by(group) %>%
mutate(ymin = min(y),
ymax = max(y),
xmin = x,
xmax = x + width / 2)
},
draw_group = function(data, panel_scales, coord) {
# Find the points for the line to go all the way around
data <- transform(data, xminv = x,
xmaxv = x + violinwidth * (xmax - x)) #利用transform函数为数据框mydata增加数据
newdata <- rbind(plyr::arrange(transform(data, x = xmaxv), -y),plyr::arrange(transform(data, x = xminv), y))
newdata_Polygon <- rbind(newdata, newdata[1,])
newdata_Polygon$colour<-NA
newdata_Path <- plyr::arrange(transform(data, x = xmaxv), -y)
ggplot2:::ggname("geom_flat_violin", grobTree(
GeomPolygon$draw_panel(newdata_Polygon, panel_scales, coord),
GeomPath$draw_panel(newdata_Path, panel_scales, coord))
)
},
draw_key = draw_key_polygon,
default_aes = aes(weight = 1, colour = "grey20", fill = "white", size = 0.5,
alpha = NA, linetype = "solid"),
required_aes = c("x", "y")
)
# "%||%" <- getFromNamespace("%||%", "ggplot2")
# "%>%" <- getFromNamespace("%>%", "magrittr")
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20）
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)+3
)
#-------------------------------------------------------------
colnames(mydata)<-c("Class", "Value")
d <- group_by(mydata, Class) %>%
summarize(mean = mean(Value),
sd = sd(Value))
ggplot(mydata, aes(Class, Value, fill=Class)) +
geom_flat_violin(position=position_nudge(x=.2)) +
geom_jitter(aes(color=Class), width=.1) +
geom_pointrange(aes(y=mean, ymin=mean-sd, ymax=mean+sd),
data=d, size=1, position=position_nudge(x=.2)) +
coord_flip() +
theme_bw() +
theme( axis.text = element_text(size=13),
axis.title = element_text(size=15),
legend.position="none")
ggplot(mydata, aes(x=Class, y=Value)) +
geom_flat_violin(aes(fill=Class),position=position_nudge(x=.25),color="black") +
geom_jitter(aes(color=Class), width=0.1) +
geom_boxplot(width=.1,position=position_nudge(x=0.25),fill="white",size=0.5) +
coord_flip() +
theme_bw() +
theme( axis.text = element_text(size=13),
axis.title = element_text(size=15),
legend.position="none")
所示。相比于
的海盗图，它显得没那么余。相比于
的小提琴图，它又省却多余的一半核密度估计曲线的同时，增加了抖动散点图。技能云雨图云雨图可以看成核密度估计曲线图、箱形图和抖动散点图的组合图表，那么就可以使用自定义的半小提琴函数geomflatviolin()、箱形图函数geom_boxplot()和抖动散点图函数geom_jitter()分别叠加实现。其中只需要将半小提琴（即核密度估计曲线）和箱形图，通过设定参数position=position_nudge(x)，将其移动到左边或上边距离X轴类别中心线的x位置，具体代码如下所示。
geom_flat_violin(aes（fill=Class),position=position_nudge(x=.25),color=“black")+geom_jitter(aes(color=Class),width=0.1)+双数据系列的箱形图、小提琴图和豆状图如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
#-------------------------------------------------(a)和(b) 多数据系列的箱型图--------------------------------------------------------
set.seed(141079)
data <- data.frame(BAI2013 = rnorm(300),
class = rep(letters[1:3], 100),
treatment = rep(c("elevated","ambient"),150))
#(a)多数据系列的箱型图
ggplot(data, aes(x = class, y = BAI2013))+
geom_boxplot(outlier.size = 1, aes(fill=factor(treatment)),
position = position_dodge(0.8),size=0.5) +
guides(fill=guide_legend(title="treatment"))+
theme_minimal()+
theme(axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="right",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#(b) 带抖动散点的多数据系列箱型图
data<-transform(data,dist_cat_n=as.numeric(class),
scat_adj=ifelse(treatment == "ambient",-0.2,0.2))
ggplot(data, aes(x =class, y = BAI2013))+
geom_boxplot(outlier.size = 0, aes(fill=factor(treatment)),
position = position_dodge(0.8),size=0.4) +
geom_jitter(aes(scat_adj+dist_cat_n, BAI2013,fill = factor(treatment)),
position=position_jitter(width=0.1,height=0),
alpha=1,
shape=21, size = 1.5)+
guides(fill=guide_legend(title="treatment"))+
theme_minimal()+
theme(axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="right",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#-------------------------------------------------(c)多数据系列的小提琴图SplitViolin------------------------------------------------
#Reference:
# https://stackoverflow.com/a/45614547
# https://gist.github.com/Karel-Kroeze/746685f5613e01ba820a31e57f87ec87
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin,
draw_group = function(self, data, ..., draw_quantiles = NULL){
# Original function by Jan Gleixner (@jan-glx)
# Adjustments by Wouter van der Bijl (@Axeman)
data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
grp <- data[1,'group']
newdata <- plyr::arrange(transform(data, x = if(grp%%2==1) xminv else xmaxv), if(grp%%2==1) y else -y)
newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
newdata[c(1,nrow(newdata)-1,nrow(newdata)), 'x'] <- round(newdata[1, 'x'])
if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <= 1))
quantiles <- create_quantile_segment_frame(data, draw_quantiles, split = TRUE, grp = grp)
aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
aesthetics$alpha <- rep(1, nrow(quantiles))
both <- cbind(quantiles, aesthetics)
quantile_grob <- GeomPath$draw_panel(both, ...)
ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
}
else {
ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
}
}
)
create_quantile_segment_frame <- function (data, draw_quantiles, split = FALSE, grp = NULL) {
dens <- cumsum(data$density)/sum(data$density)
ecdf <- stats::approxfun(dens, data$y)
ys <- ecdf(draw_quantiles)
violin.xminvs <- (stats::approxfun(data$y, data$xminv))(ys)
violin.xmaxvs <- (stats::approxfun(data$y, data$xmaxv))(ys)
violin.xs <- (stats::approxfun(data$y, data$x))(ys)
if (grp %% 2 == 0) {
data.frame(x = ggplot2:::interleave(violin.xs, violin.xmaxvs),
y = rep(ys, each = 2), group = rep(ys, each = 2))
} else {
data.frame(x = ggplot2:::interleave(violin.xminvs, violin.xs),
y = rep(ys, each = 2), group = rep(ys, each = 2))
}
}
geom_split_violin <- function (mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, show.legend = NA, inherit.aes = TRUE) {
layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, position = position, show.legend = show.legend, inherit.aes = inherit.aes, params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}
data<-transform(data,dist_cat_n=as.numeric(class),
scat_adj=ifelse(treatment == "ambient",-0.15,0.15))
#data$scat_adj[data$treatment == "ambient"] <- -0.15
#data$scat_adj[data$treatment == "elevated"] <- 0.15
ggplot(data, aes(x = class, y = BAI2013,fill=factor(treatment)))+
geom_split_violin(draw_quantiles = 0.5,trim = FALSE)+
geom_jitter(aes(scat_adj+dist_cat_n, BAI2013,fill = factor(treatment)),
position=position_jitter(width=0.1,height=0),
alpha=1,
shape=21, size = 1)+
guides(fill=guide_legend(title="treatment"))+
theme_minimal()+
theme(axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="right",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#----------------------------------------------------(d) 多数据系列的豆状图-----------------------------
library(beanplot)
par(mai=c(0.5,0.5,0.25,1.2))
beanplot(BAI2013 ~treatment*class, data,col = list("#FF6B5E", "#00C3C2"),
side = "both",xlab ="Class",ylab ="value")
legend(x=3.7,y=1.5 ,xpd=TRUE,bty="n",c("ambient", "elevated"),
fill = c("#FF6B5E", "#00C3C2"),title="treatment")
所示。双数据系列的箱形图可以使用geom_boxplot()函数，只需要将两组的变量映射到箱形的填充颜色（fill），另外可以使用position=position_dodge（width）控制箱形之间的间隔，如
所示。在
基础上，可以再使用geom_jitter()函数添加抖动散点图，可以通过position=position_jiter（width，height语句使散点沿着箱形图的中心线分布，如
所示。
是双数据系列的小提琴图，它并非像双数据系列的箱形图一样，同一个类别下，两个小提琴图。这是因为小提琴图本身就是由两个左右对称的核密度估计曲线图构成的。所以对于双数据系列小提琴图，我们只需要保留两个小提琴图的各一半，使左边为一个数据的核密度估计曲线图，右边为另一个数据的核密度估计曲线图。
由于ggplot2包并未提供这样的函数，所以我们可以通过自定义双数据系列小提琴图的绘制函数geom_split_violin()实现。在此基础上，再使用geom_jitter()函数添加抖动散点图。
是双数据系列的豆状图，可以使用beanplot包的beanplot()函数就可以直接实现。跟
（c)小提琴图表达的数据信息基本一致。中中(b）带抖动散点的双数据系列箱形图2.50.0（）带抖动散点的双数据系列小提琴图技能双数据系列箱形图双数据系列箱形图的核心代码如下所示
双数据系列的箱形图geom_boxplot（outlier.size=1,aes（fill=factor（treatment),position =position_dodge(0.8),size=0.5)+#
带抖动散点的多数据系列箱形图data<-transform(data,dist_cat_n=as.numeric(class),scat_adj=ifelse(treatment =="ambient",-0.2,0.2)geom_boxplot（outlier.size=0.aes（fill=factor（treatment),position=position_dodge(0.8),size=0.4)+geom_jitter(aes（scat_adj+dist_cat_n,BAl2013,fill =factor（treatment)),position=position_jitter（width=0.1,height=0)
为带抖动散点的双数据系列小提琴图，需要使用自定义的函数geom_split_violin()实现，它可以将两个小提琴图各取一半，并拼接在一起，具体实现代码如下所示
）。它与
所示组合图表类似，但不完全相同。子母图是在主图的基础上，再添加子图。4.54.03.53.02.52.0（a)二维散点图与箱形图的组合图表技能子母图ggplot2包也可以实现绘制子母图，通过viewport()函数来实现，viewport()是grid绘图体系用于排版的函数（ggplot2包是基于grid绘图原理设计的）：viewport()函数主要的参数有4个，x和y设置中心位点相对于父图层的位置，width和height设置子图形的大小，如
所示。
所示的子母图的具体实现代码如下所示。宽度高度p1<-ggplot(iris,aes(Sepal.Length,Sepal.Width,fill=Species))+geom_point（size=4,shape=21,color=“black")+p2<-ggplot(iris,aes（Species,Sepal.Width,fill =Species)+subvp<-viewport（x=0.78.y=0.38,width=0.4.height=0.5)有时候，我们需要对箱形图进一步添加显著性标签，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,3,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(3, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20
k <- rJohnson(100, findParams(3, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#------------------------------------------图5-2-11 带显著性标签的箱型图(a)-----------------------------------------
library(ggpubr)
palette<-c(brewer.pal(7,"Set2")[c(1,2,4,5)])
ggboxplot(mydata, x = "Class", y = "Value",
fill = "Class", palette = palette,
add = "none",size=0.5,add.params = list(size = 0.25))+
geom_hline(yintercept = mean(mydata$Value), linetype = 2)+ #添加均值线
stat_compare_means(method = "anova", label.x=0.8,label.y = 7.8)+ # 添加全部数据的annova 方法的p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.", hide.ns = TRUE,label.y = 8) + # 添加每组变量与全部数据的显著性
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------------------Method1:图5-2-11 带显著性标签的箱型图(b)-----------------------------------------
compaired <- list(c("n", "s"),
c("n","k"),
c("n","mm"),
c("s","k"))
ggboxplot(mydata, x = "Class", y = "Value",
fill = "Class", palette = palette,
add = "jitter",size=0.5)+
stat_compare_means(comparisons = compaired,method = "wilcox.test")+ # 添加每两组变量的显著性
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------------------Method2:-图5-2-11 带显著性标签的箱型图(b)-----------------------------------------
ggplot(mydata, aes(Class, Value))+
geom_boxplot(aes(fill = Class),notch = FALSE,outlier.alpha =1) +
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
geom_signif(comparisons = compaired,
step_increase = 0.1,
map_signif_level = F,
test = wilcox.test)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
所示。R中常用的比较方法主要如表5-2-1所示。表5-2-1R中常用的比较方法方法R函数描述T检验，比较两组（参数）Wilcoxon符号秩检验，比较两组（非参数）aov()或anova()方差检验，比较多组（参数）Kruskal-Wallis检验，比较多组（非参数）10.00.330.230.000487.52.50.0技能带显著性标签的箱形图ggpubr包中的两个函数：compare_means()可以进行一组或多组间的比较；stat_compare_mean()自动添加p-value、显著性标签到ggplot2绘制的图表。
带显著性标签的箱形图的具体代码如下所示。ggpubr包用于出版物图表的绘制。HadleyWickham创建的可视化包ggplot2可以流畅地进行优美的可视化，但是如果要通过ggplot2定制一套图形，尤其是适用于杂志期刊等出版物的图形，对于那些没有深入了解ggplot2的人来说就有点困难了，ggplot2的部分语法是很晦涩的。
为此AlboukadelKassambara创建了基于ggplot2的可视化包ggpubr，用于绘制符合出版物要求的图形
带显著性标签的箱形图add=“none"size=0.5,add.params=list（size=0.25))+#添加均值线#添加全部数据的annova方法的p-valuestat_compare_means（label =“p.signif",method =“t.test",ref.group ="all",hide.ns =TRUE,label.y= 8) +#添加每组变量与全部数据的显著性#
带显著性标签的箱形图compaired<-list（c("n","s"),c("n","k"),c("n""mm"),c("s","k")stat_compare_means(comparisons=compaired,method=“wilcox.test")有时候，我们还需要比较两个成对样本（pairedsample），比如在对高血压的研究中，在研究开始会测量所有病人的血压，在治疗之后再次测量血压。这样，每个主体有两个测量值，它们通常称为之前测量值和之后测量值，这就是成对样本。成对样本T检验一般是比较单独一组的两个变量的平均值。
此过程计算每个个案的两个变量的值之间的差值，并检验平均差值是否非0（见
#EasyCharts团队出品，如有商用必究，
#如需使用与深入学习，请联系微信：EasyCharts
library(ggplot2)
#---------------------------------------------Method1：图5-2-14 带连接线的双箱型图--------------------------------------------------------
library(ggpubr)
set.seed(141079)
data <- data.frame(BAI2013 = rnorm(60),
class = rep(rep(letters[1:3], each=10),2),
treatment = rep(c("elevated","ambient"),each=30),
index=rep(seq(1,30),2))
palette<-c(brewer.pal(7,"Set2")[c(1,2,4,5)])
ggpaired(data, x = "treatment", y = "BAI2013",
fill = "treatment", palette = palette,
line.color = "grey50", line.size = 0.15, point.size = 1.5,width=0.6,
facet.by = "class", short.panel.labs = FALSE)+
stat_compare_means(paired = TRUE)+
theme_minimal()+
theme(strip.background = element_rect(fill="grey90"),
strip.text = element_text(size=13,face="plain",color="black"),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid=element_blank(),
legend.position="none",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#-------------------------------------------Method2：图5-2-14 带连接线的双箱型图--------------------------------------------------
library(RColorBrewer)
library(reshape2)
library(ggforce)
library(dplyr)
set.seed(141079)
df_point <- data.frame(BAI2013 = rnorm(60),
class = rep(rep(letters[1:3], each=10),2),
treatment = rep(c("elevated","ambient"),each=30),
index=rep(seq(1,30),2))
type<-as.character(unique(df_point$class))
df_bezier<-data.frame(matrix(ncol = 4, nrow = 0))
colnames(df_bezier)<-c("index","treatment","class","value")
for (i in 1:length(type)){
data0<-df_point[df_point$class==type[i],]
data1<-split(data0,data0$treatment)
data2<-data.frame(ambient=data1$ambient[,1],
elevated=data1$elevated[,1],
index=data1$ambient[,4])
colnames(data2)<-c(1,2,"index")
data2$'1.3'<-data2$'1'
data2$'1.7'<-data2$'2'
data3<-melt(data2,id="index",variable.name ="treatment")
data3$treatment<-as.numeric((as.character(data3$treatment)))
data4<-arrange(data3,index,treatment)
data4$class<-type[i]
df_bezier<-rbind(df_bezier,data4)
}
ggplot()+
geom_boxplot(data=df_point,aes(x = factor(treatment), y = BAI2013,fill=factor(treatment)),
width=0.35,position = position_dodge(0),size=0.5,outlier.size = 0) +
geom_point(data=df_point,aes(x = factor(treatment), y = BAI2013,fill=factor(treatment)),
shape=21,colour="black",size=2)+
#geom_line(data=df_point,aes(x = factor(treatment), y = BAI2013,group=index),
# size=0.25,colour="grey20")+
geom_bezier(data=df_bezier,aes(x= treatment, y = value, group = index,linetype = 'cubic'),
size=0.25,colour="grey20") +
scale_fill_manual(values=brewer.pal(7,"Set2")[c(5,2)])+
facet_grid(.~class)+
labs(x="treatment",y="Value")+
theme_minimal()+
theme(strip.background = element_rect(fill="grey90"),
strip.text = element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="none",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
）。图：5-2-14带连接线的双箱形图技能带连接线的双箱形图对于成对样本，我们可以使用ggpubr包的ggpaired()函数实现可视化，但是每对样本数据之间是使用直线连接的，这就导致数据可视化效果并不美观，所以我们可以自定义绘图，使用ggforce包的geom_bezier()函数，从而用光滑的贝塞尔曲线连接两点，如
所示。其关键在于构造贝塞尔连接曲线的数据集，具体代码如下所示
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(ggstance) #devtools::install_github("lionel-/ggstance")
color<-brewer.pal(7,"Set2")[c(1,2,4,5)]
set.seed(141079)
data <- data.frame(BAI2013 = rnorm(300), class = rep(letters[1:3], 100),
treatment = rep(c("elevated","ambient"),150))
data<-transform(data,dist_cat_n=as.numeric(class), scat_adj=ifelse(treatment == "ambient",-0.2,0.2))
#--------------------------------------------------图5-2-13 水平显示的箱型图(a)-----------------------------------------------
ggplot(data, aes(class,BAI2013))+
geom_boxplot(aes(fill=factor(treatment)),
size=0.5,outlier.size = 1,
position = position_dodge(0.8)) +
guides(fill=guide_legend(title="treatment"))+
theme_minimal()+
coord_flip()+
theme(axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="right",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
#----------------------------------------------------图5-2-13 水平显示的箱型图(b)-----------------------------------------------
ggplot(data, aes(BAI2013,class))+
geom_boxploth(aes(fill=factor(treatment)),
size=0.5,outlier.size = 1,
position =position_dodgev(0.8)) +
guides(fill=guide_legend(title="treatment"))+
theme_minimal()+
theme(axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=11,face="plain",color="black"),
panel.background=element_rect(colour="black",fill=NA),
panel.grid.minor=element_blank(),
legend.position="right",
legend.background=element_rect(colour=NA,fill=NA),
axis.ticks=element_line(colour="black"))
所示。虽然箱形图部分实现了水平翻转，但是右边的图例（legend）部分还是竖直的。这时，我们只需要把ggplot2包的geom_box()函数替换成ggstance包的geom_boxploth()函数，就可以实现
所示的效果。箱形图的中值排序显示：排序显示数据对更快地发现数据规律和获取数据信息尤为重要。对应X轴为类别向量时，最好将箱形图按中值降序后显示，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
library(SuppDists) #提供rJohnson()函数
set.seed(141079)
# Generate sample data -------------------------------------------------------
#findParams函数参考：https://github.com/hadley/boxplots-paper
findParams <- function(mu, sigma, skew, kurt) {
value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma),
as.double(skew), as.double(kurt - 3), gamma = double(1),
delta = double(1), xi = double(1), lambda = double(1),
type = integer(1), PACKAGE = "SuppDists")
list(gamma = value$gamma, delta = value$delta,
xi = value$xi, lambda = value$lambda,
type = c("SN", "SL", "SU", "SB")[value$type])
}
# 均值为3，标准差为1的正态分布
n <- rnorm(100,8,1)
# Johnson分布的偏斜度2.2和峰度13
s <- rJohnson(100, findParams(4, 1, 2., 13.1))
# Johnson分布的偏斜度0和峰度20
k <- rJohnson(100, findParams(10, 1, 2.2, 20))
# 两个峰的均值μ1，μ2分别为1.89和3.79，σ1 = σ2 =0.31
mm <- rnorm(100, rep(c(2, 4), each = 50) * sqrt(0.9), sqrt(0.1))
mydata <- data.frame(
Class = factor(rep(c("n", "s", "k", "mm"), each = 100),
c("n", "s", "k", "mm")),
Value = c(n, s, k, mm)
)
#write.csv(mydata,'Boxplot_Sort_Data.csv')
#--------------------------------未排序----------------------------------------------------
ggplot(mydata, aes(Class,Value))+
geom_boxplot(aes(fill = Class),notch = FALSE,outlier.alpha =1) +
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
scale_y_continuous(breaks=seq(0,15,3))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
#------------------------------降序序处理-------------------------------------------------------
Order_Class<-with(mydata,reorder(Class,Value,median))
Order_Class<-factor(Order_Class,levels=rev(levels(Order_Class)))
ggplot(mydata, aes(Order_Class,Value))+
geom_boxplot(aes(fill = Class),notch = FALSE,outlier.alpha =1) +
scale_fill_manual(values=c(brewer.pal(7,"Set2")[c(1,2,4,5)]))+
scale_y_continuous(breaks=seq(0,15,3))+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="none"
)
所示。技能中值排序显示的箱形图先使用reorder()函数根据中值（median）对数据框mydata排序，然后改变因子向量的顺序，使因子向量的水平（level）按其中值的降序排列，最后使用ggplot2包的geomboxplot()函数绘制即可，具体代码如下所示
5.3 ：二维统计直方图和二维核密度估计图
一句话：两个变量的联合分布：二维直方图和二维核密度。
5.3.1 二维统计直方图
一句话：二维直方图：X/Y都分组计数，看二维数据聚集区域。
二维统计直方图主要针对二维数据的统计分析，X-Y轴变量为数值型。首先要从X轴和Y轴变量数据分别找出它的最大值和最小值，然后确定一个区间，使其包含全部测量数据，将区间分成若干小区间[x：X+w，Y：Y+w]（其中，w为最小区间的大小，（x，Y）为第n个区间的始点），统计测量结果出现在各小区间的频数M。在平面直角坐标系中，X轴和Y轴分别标出每个组的端点，每个方块（bin）的颜色代表对应的频数，一般我们也称这样的统计图为二维频数分布直方图（见
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(RColorBrewer)
colormap<- rev(brewer.pal(11,'Spectral'))
# Create normally distributed data for plotting
x1 <- rnorm(mean=1.5, 5000)
y1 <- rnorm(mean=1.6, 5000)
x2 <- rnorm(mean=2.5, 5000)
y2 <- rnorm(mean=2.2, 5000)
x<-c(x1,x2)
y<-c(y1,y2)
df <- data.frame(x,y)
#--------------------------图5-3-1 不同类型的二维统计直方图------------------
ggplot(df, aes(x,y))+
stat_bin2d(bins=40) + scale_fill_gradientn(colours=colormap)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="right"
)
ggplot(df, aes(x,y))+
stat_binhex(bins=40) + scale_fill_gradientn(colours=colormap)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="right"
)
5.3.2 二维核密度估计图
一句话：二维核密度估计图：二维平滑密度，比二维直方图更顺滑。
关于核密度估计，我们在5.1.2节有过介绍，此处不再赘述。二维核密度估计图如
#EasyChartsŶӳƷñؾ
#ʹѧϰϵţEasyCharts
library(ggplot2)
library(RColorBrewer)
colormap<- rev(brewer.pal(11,'Spectral'))
# Create normally distributed data for plotting
x1 <- rnorm(mean=1.5, 5000)
y1 <- rnorm(mean=1.6, 5000)
x2 <- rnorm(mean=2.5, 5000)
y2 <- rnorm(mean=2.2, 5000)
x<-c(x1,x2)
y<-c(y1,y2)
df <- data.frame(x,y)
#------------------------------------ͼ5-3-2 ͬ͵Ķάܶͳͼ-----------------
ggplot(df, aes(x,y))+
stat_density_2d(geom ="raster",aes(fill = ..density..),contour = F)+# "polygon")+#geom_raster(aes(fill = density)) +
scale_fill_gradientn(colours=colormap)+#, trans="log"scale_fill_gradientn(colours=c("#CEF5FF","#00B8E5","#005C72"),name = "Frequency",na.value=NA)+
#scale_fill_gradientn(colours=c(brewer.pal(7,"Set2")[3],"white",brewer.pal(7,"Set2")[2]),na.value=NA)+
#geom_contour(acolour = "white") +
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="right"
)
#---------------------------
ggplot(df, aes(x, y)) +
stat_density2d(geom ="polygon",aes(fill = ..level..),bins=30 )+#alpha=..level..,aes( fill=..level..), size=2, bins=10, geom="polygon") +
#stat_density_2d(geom = "point", aes(size = ..density..), n = 20, contour = FALSE)
scale_fill_gradientn(colours=colormap)+#scale_fill_gradient(low = "yellow", high = "red") +
#scale_alpha(range = c(0.00, 0.5), guide = FALSE) +
#geom_density2d( colour=NA,bins=30) +#
#geom_point() +
guides(alpha=FALSE) +
xlim(-2,6)+
ylim(-2,6)+
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position="right"
)
所示。0.1250.1250.1000.1000.0750.075技能二维统计直方图和二维核密度估计图对于二维统计直方图，R中ggplot2包的geom_bin2d0函数和geomhex()函数分别可以绘制
，参数bins为在X轴和Y轴变量分别设定的区间数目。对于二维核密度估计图，R中ggplot2包的stat_density_2d0函数可以绘制，其中geom=raster”或者polygon分别对应
，具体代码如下所示
（b)geom_bin2d（bins=40.na.rm=TRUE)+#对应
（a)scale_fill_gradientn(colours=colormap)+stat_density_2d（geom=raster"aes（fill=density.）,contour=F)+#对应
（a)#stat_density_2d（geom="polygon",aes（fill=.level.),bins=30)+#对应
scale_fill_gradientn(colours= colormap)+技能三维统计分布图对于
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(plot3D)
N<-300
x1 <- rnorm(mean=1.5, N)
y1 <- rnorm(mean=1.6, N)
x2 <- rnorm(mean=2.5, N)
y2 <- rnorm(mean=2.2, N)
data <- data.frame(x=c(x1,x2),y=c(y1,y2))
#图5-3-3 (a) 三维统计直方图
library(gplots) #提供hist2d()函数
df_hist<-hist2d(df$x,df$y, nbins=30)
pmar <- par(mar = c(5.1, 4.1, 4.1, 6.1))
hist3D(x=df_hist$x,y=df_hist$y,z=df_hist$counts,
col = colormap, border = "black",space=0,alpha = 1,lwd=0.1,
xlab = "x", ylab = "y",zlab = "Count", clab="Count",
ticktype = "detailed",bty = "f",box = TRUE,#cex.axis= 1e-09,
theta = 65, phi = 20, d=3,
colkey = list(length = 0.5, width = 1))
#图5-3-3 (b) 三维核密度估计图
library(MASS) #提供kde2d ()函数
df_density <- kde2d(df$x,df$y, n = 50, h = c(width.SJ(df$x), width.SJ(df$y)))
pmar <- par(mar = c(5.1, 4.1, 4.1, 6.1))
persp3D (df_density$x, df_density$y, df_density$z,
theta = 60, phi = 20, d=3,
col = colormap, border = "black", lwd=0.1,
bty = "f",box = TRUE,ticktype = "detailed",
xlab = "x", ylab = "y",zlab = "desnity",clab="desnity",
colkey = list(length = 0.5, width = 1))
所示的三维统计直方图，其绘制方法是先使用gplots包的hist2d0函数求二维统计直方图数值，其中bins表示X轴和Y轴方向的箱形总数，最后使用plot3D包的hist3D0函数绘制三维柱形图。对于
所示的三维核密度估计图，其绘制方法是先使用MASS包的kde2d0函数计算二维核密度估计，其中h为X和Y轴方向的带宽，最后使用plot3D包的persp3DO函数绘制三维曲面图。0.100.120.100.080.060.040.000.020.00
所示图表的实现代码如下所示。#
三维统计直方图df_hist<-hist2d(df$x,df$y.nbins=30)pmar<-par（mar=c(5.1,4.1,4.1,6.1)hist3D(x=df_hist$x.y=df_hist$y.z=df_hist$counts,col =colormap,border=“black"space=0,alpha =1,lwd=0.1,xlab="x".ylab="y”zlab=“Count".clab="Count”,ticktype ="detailed"bty="f"box=TRUE,#cex.axis=1e-09colkey=list（length =0.5.width =1))#
三维核密度估计图pmar<-par（mar=c(5.1,4.1,4.1,6.1))persp3D（df_density$x,df_density$y.df_density$zcol=colormap,border=“black",Iwd=0.1,bty=“f",box=TRUE,ticktype=“detailed”xlab=“x",ylab=“y”zlab=“desnity"clab="desnity”colkey=list（length =0.5,width =1))二维与一维统计分布组合图：我们还可以将二维统计直方图和二维核密度估计图，结合一维的统计分布图表一起展示，更加详细地揭示数据的分布情况，如
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
library(ellipse)
library(gridExtra)
library(plyr)
library(RColorBrewer)
Colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
N<-300
x1 <- rnorm(mean=1.5, N)
y1 <- rnorm(mean=1.6, N)
x2 <- rnorm(mean=2.5, N)
y2 <- rnorm(mean=2.2, N)
data <- data.frame(x=c(x1,x2),y=c(y1,y2))
#-------------------------二维直方图+一维直方图------------------------------------------------------------
# 绘制上边的直方图，并将各种标注去除
hist_top <- ggplot()+
geom_histogram(aes(data$x),colour='black',fill='#5E4FA2',binwidth = 0.3)+
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.line=element_blank())
# 同样绘制右边的直方图
hist_right <- ggplot()+
geom_histogram(aes(data$y),colour='black',fill='#5E4FA2',binwidth = 0.3)+
theme(panel.background=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.line=element_blank())+
coord_flip()
#ggplot(diamonds, aes(carat, price))
scatter<-ggplot(data, aes(x,y)) +
stat_binhex(bins = 15,na.rm=TRUE,color="black")+#colour="black",
scale_fill_gradientn(colours=Colormap)+#, trans="log"
#geom_point(colour="white",size=1,shape=21) +
theme_classic()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position=c(0.10,0.80),
legend.background=element_blank()
)
# 最终的组合
grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(4,1), heights=c(1,4))
#----------------------------二维核密度估计图+一维核密度估计图--------------------------
# 绘制上边的直方图，并将各种标注去除
hist_top <- ggplot(data, aes(x)) +
geom_density(colour="black",fill='#5E4FA2',size=0.25)+
theme_void()
# 同样绘制右边的直方图
hist_right <- ggplot(data, aes(y)) +
geom_density(colour="black",fill='#5E4FA2',size=0.25)+
theme_void()+
coord_flip()
scatter<-ggplot(data, aes(x, y)) +
stat_density2d(geom ="polygon",aes(fill = ..level..),bins=30 )+#alpha=..level..,aes( fill=..level..), size=2, bins=10, geom="polygon") +
scale_fill_gradientn(colours=Colormap)+#, trans="log"
#geom_point(size=1) +
theme_minimal()+
theme(panel.background=element_rect(fill="white",colour="black",size=0.25),
#panel.grid.major = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#panel.grid.minor = element_line(colour = "grey60",size=.25,linetype ="dotted" ),
#text=element_text(size=15),
#plot.title=element_text(size=15,family="myfont",hjust=.5),
axis.line=element_line(colour="black",size=0.25),
axis.title=element_text(size=13,face="plain",color="black"),
axis.text = element_text(size=12,face="plain",color="black"),
legend.position=c(0.9,0.22),
legend.background=element_blank()
)
# 最终的组合
grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(4,1), heights=c(1,4))
所示。0.100.05技能二维与一维统计分布组合图R中gridExtra包的grid.arrange()函数可以实现ggplot2包绘制的一维和二维统计分布图的组合，具体实现代码如下所示
5.4 金字塔图和镜面图
一句话：金字塔图：背靠背直方图，人口结构等对比型分布。
金字塔图通常用来理解人口结构，也被称为人口金字塔（populationpyramid）图，是彼此背靠背的一对直方图，通常用于显示所有年龄组和男女人口的分布情况。X轴表示人口数量，Y轴列出年龄组别。人口金字塔图最适合用来检测人口模式的变化或差异。
多个人口金字塔图放在一起可用于比较各国或不同群体之间的人口模式。举个例子，底部较宽、顶部狭窄的人口金字塔图表示该群体具有很高的生育率和死亡率：相反，顶部较宽、底部狭窄的人口金字塔图代表出现人口老龄化，而且生育率低。除此之外，人口金字塔图也可用来推测人口的未来发展。
如果人口出现老龄化，而且生育率低，最终会导致因没有足够后代照顾老人的社会问题。其他理论包括“青年膨胀”，即若社会存在大量16~30岁的青年（特别是男性），则容易导致社会动荡、战争和恐怖主义。因此，人口金字塔图对生态学、社会学和经济学等领域都相当有用。
技能金字塔图
#EasyCharts团队出品，
#如有问题修正与深入学习，可联系微信：EasyCharts
library(ggplot2)
#-----------------------------------------------(a1)-------------------------------------------
df<-read.csv("Population_Pyramid_Data.csv",header=TRUE)
df[df$gender == "female",]$pop<--df[df$gender == "female",]$pop
df$age<-factor(df$age,levels=df$age[seq(1,nrow(df)/2,1)])
ggplot(data = df, aes(x =age , y = pop, fill = gender)) +
geom_bar(stat = "identity",position = "identity",color="black",size=0.25) +
scale_y_continuous(labels = abs, limits = c(-400, 400), breaks = seq(-400, 400, 100)) +
coord_flip() +
theme_light()+
theme(
#axis.text.x = element_text(angle=60, hjust=1),
panel.grid.minor=element_blank(),
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=10,face="plain",color="black"),
legend.title=element_text(size=14,face="plain",color="black"),
legend.text=element_text(size=12,face="plain",color="black"),
legend.background=element_blank(),
legend.position = c(0.9,0.88)
)
#------------------------------------(a2)------------------------------------
ggplot(data = df, aes(x =age , y = pop, fill = gender)) +
geom_bar(stat = "identity",position = "identity",color="black",size=0.25) +
scale_y_continuous(labels = abs, limits = c(-400, 400), breaks = seq(-400, 400, 100)) +
theme_light()+
theme(
axis.text.x = element_text(angle=60, hjust=1),
panel.grid.minor=element_blank(),
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=10,face="plain",color="black"),
legend.title=element_text(size=14,face="plain",color="black"),
legend.text=element_text(size=12,face="plain",color="black"),
legend.background=element_blank(),
legend.position = c(0.9,0.88)
)
#-------------------------------------------(b1)--------------------------------------------------------
#df<-read.csv("Population_Pyramid_Data.csv",header=TRUE)
#df[df$gender == "female",]$pop<--df[df$gender == "female",]$pop
df$age_x<-rep(seq(0, 100,5),2)
ggplot(data = df, aes(x =age_x , y = pop, fill = gender)) +
geom_area(stat = "identity", position = "identity",color="black",size=0.25) +
scale_fill_manual(values=c("#36BED9","#FBAD01"))+
coord_flip() +
scale_y_continuous(labels = abs, limits = c(-400, 400), breaks = seq(-400, 400, 100)) +
scale_x_continuous(breaks = seq(0, 100, 5),labels=df$age[seq(1,nrow(df)/2,1)])+
theme_light()+
theme(
panel.grid.minor=element_blank(),
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=10,face="plain",color="black"),
legend.title=element_text(size=14,face="plain",color="black"),
legend.text=element_text(size=12,face="plain",color="black"),
legend.background=element_blank(),
legend.position = c(0.9,0.88)
)
#-------------------------------------------(b2)------------------------------------------------------------------
ggplot(data = df, aes(x =age_x , y = pop, fill = gender)) +
geom_area(stat = "identity", position = "identity",color="black",size=0.25) +
scale_fill_manual(values=c("#36BED9","#FBAD01"))+
#coord_flip() +
scale_y_continuous(labels = abs, limits = c(-400, 400), breaks = seq(-400, 400, 100)) +
scale_x_continuous(breaks = seq(0, 100, 5),labels=df$age[seq(1,nrow(df)/2,1)])+
theme_light()+
theme(
axis.text.x = element_text(angle=60, hjust=1),
panel.grid.minor=element_blank(),
#text=element_text(size=15,face="plain",color="black"),
axis.title=element_text(size=15,face="plain",color="black"),
axis.text = element_text(size=10,face="plain",color="black"),
legend.title=element_text(size=14,face="plain",color="black"),
legend.text=element_text(size=12,face="plain",color="black"),
legend.background=element_blank(),
legend.position = c(0.9,0.88)
)
（a2）其实就是由两个不同数据系列的柱形图组成的，R中ggplot2包的同数据系列的面积图组成的，R中ggplot2包的geomarea()函数可以绘制基于面积的镜面图和金字塔图，只是由于面积图的X轴只能为数值型，所以需要构造X轴辅助数据，然后将X轴坐标标签替换成年龄分段。
（b1）所示图表的实现代码如下所示
（a1）直方图类型的金字塔图geom_bar(stat=“identity"position=“identity"color="black"size=0.25)+scale_y_continuous（labels=abs,limits=c(-400,400),breaks=seq(-400,400,100)+#
（b1）面积图类型的金字塔图df$age_x<-rep(seq（0,100,5).2)geom_area(stat=“identity".position=“identity"color=“black"size=0.25)+scale_fill_manual(values=c(#36BED9""#FBAD01"))+scale_x_continuous（breaks=seq（0,100,5),labels=df$age[seq（1,nrow（df)/2,1)])+