R 案例｜绘制不同分布的 QQ 图
R 案例｜绘制不同分布的 QQ 图
简介
论文中需要绘制数据对于不同分布假定下的 QQ 图。这里小编主要是使用 qqplotr 包进行绘制，参考的博客：An Introduction to qqplotr[1]。
简单版本
绘制正态分布的 QQ 图
对于经典的正态分布的 QQ 图，大家可能并不陌生，并且在网上可以找到很多“搬运”的中文推文。但是解释的都不是很清楚。
这里我以这篇博客中的某个例子为例，进行介绍：
1. 加载包
library(qqplotr)
library(ggplot2)
2. 随机产生数据
set.seed(0)
smp <- data.frame(norm = rnorm(100))
3. 绘制正态分布的 QQ 图
gg <- ggplot(data = smp, mapping = aes(sample = norm)) +
stat_qq_band() +
stat_qq_line() +
stat_qq_point() +
labs(x = "Theoretical Quantiles", y = "Sample Quantiles")
gg
- 拓展
这里做一个简单拓展，如果你想使用不同的置信带构造置信区间，可以用参数 bandType 进行调整。下面代码给出三种不同方法构造置信区间的结果。并且使用 viridis 包，对其进行配色修改。
library(viridis)
gg <- ggplot(data = smp, mapping = aes(sample = norm)) +
geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.9) +
geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.9) +
geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.9) +
geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.9) +
stat_qq_line() +
stat_qq_point() +
labs(x = "Theoretical Quantiles", y = "Sample Quantiles") +
scale_fill_viridis(discrete = T,direction = -1)
gg
进阶版本
读者绘制正态分布的 QQ 图，还是比较简单。但是如果是其他分布的情况呢？
这里以一个可靠性数据为例子，该数据来源于文献：Badar, M. G., Priest, A. M. (1982). Statistical aspects of fiber and bundle strength in hybrid composites. In: Hayashi, T., Kawata, K., Umekawa, S., eds. Progress in Science and Engineering Composites. Tokyo: ICCM-IV, pp. 1129–1136。
data = data.frame('y' = c(1.339, 1.434, 1.549, 1.574 ,1.589, 1.613, 1.746 ,1.753, 1.764 ,1.807, 1.812, 1.84, 1.852, 1.852, 1.862, 1.864, 1.931, 1.952, 1.974, 2.019, 2.051, 2.055, 2.058 ,2.088, 2.125, 2.162, 2.171, 2.172 ,2.18 ,2.194 ,2.211 ,2.27, 2.272, 2.28, 2.299, 2.308, 2.335 ,2.349 ,2.356 ,2.386, 2.39, 2.41, 2.43, 2.431, 2.458, 2.471, 2.497, 2.514 ,2.558, 2.577, 2.593, 2.601, 2.604, 2.62 ,2.633, 2.67, 2.682, 2.699, 2.705, 2.735, 2.785, 2.785,3.02, 3.042, 3.116, 3.174))
绘制指数分布的 QQ 图
这里先绘制其指数分布的 QQ 图。根据指数函数参数拟合该数据之后，得到 rate =2.2867，并将其保存到 list 中。
具体如何拟合，读者自行搜索 R 包中的相关函数。
其他代码基本不变，主要是将 stat_qq_line() 和 stat_qq_point() 中的分布设定下，参数设定下。
# exponential distribution
dp <- list(rate = 2.2867)
di <- "exp"
p1 = ggplot(data = data, mapping = aes(sample = y)) +
stat_qq_band(distribution = di, dparams = dp) +
stat_qq_line(distribution = di, dparams = dp) +
stat_qq_point(distribution = di, dparams = dp) +
labs(x = "Theoretical Quantiles", y = "Sample Quantiles") +
theme_bw() +
theme(panel.grid = element_blank())
p1
绘制威布尔分布的 QQ 图
同理，将该数据应用到威布尔分布中。结果如下：
# weibull distribution
di <- "weibull" # exponential distribution
dp <- list(shape=5.4766,scale=2.4113)
p2 = ggplot(data = data, mapping = aes(sample = y)) +
stat_qq_band(distribution = di, dparams = dp) +
stat_qq_line(distribution = di, dparams = dp) +
stat_qq_point(distribution = di, dparams = dp) +
labs(x = "Theoretical Quantiles", y = "Sample Quantiles") +
theme_bw() +
theme(panel.grid = element_blank())
p2
可以看到该数据集更适合使用 Weibull 分布进行拟合，指数分布拟合的结果较差些。读者可以使用其他分布进行拟合，并比较对应的 QQ 图，寻找最合适的分布。
然后把这些 QQ 图 合并到一起，通过可视化直观的进行比较。
这里使用 cowplot[2] 包，将两图进行合并。小编对该包的介绍做过几期，可见：cowplot包：用R添加水印。其他合并的方式还有：R可视乎｜合并多幅图形。
# 合并两幅 QQ 图
library(cowplot)
plot_grid(p1, p2, ncol = 2, nrow = 1)
参考资料
[1]
An Introduction to qqplotr: https://cran.r-project.org/web/packages/qqplotr/vignettes/introduction.html
[2]
cowplot: https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
腾讯云开发者
