简介
最近科研论文的审稿意见，需要对数据的拟合情况进行说明。分位数-分位数图[1]（Quantile-Quantile plot，简称 QQ 图）是一种不错的选择。
研究动机
在 R 语言中常用 qqnorm() 绘制 QQ 图的散点图，qqline() 用于在图上添加一条参考直线。可支持与一些常见的理论分布进行比较，例如：正态分布 ("norm")、指数分布 ("exp")、伽马分布 ("gamma")、威布尔分布 ("weibull")、对数正态分布 ("lnorm")等。
但是绘制出来的图形可能不是很美观，如下所示：
本文贡献点
当然小编以前使用 ggplot 绘制过：绘制不同分布的 QQ 图，但内容不包括逆高斯分布以及其他比较“新”的分布。
基于这个痛点，本文小编将根据 QQ 图的定义，通过 ggplot2[2] 复现自己常用的三个分布（正态、伽马和逆高斯分布）的 QQ 图。
题外话：读者可以根据本文思路，拓展到任何自己感兴趣的分布。
教程
理论介绍
介绍：一种用于检查数据是否符合某个理论分布的统计图。
基本思想：比较观察数据的分位数与理论分布的分位数。
判断依据：如果数据符合理论分布，那么 QQ 图上的点将近似落在一条直线上。
绘制步骤：
- 排序数据：将观察数据按照大小进行排序。
- 计算分位数：计算每个观察数据点对应的分位数。这可以通过排名除以总数据点的数量来实现。
- 计算理论分位数：使用所选择的理论分布的累积分布函数，计算对应于每个观察数据点的理论分位数。
- 绘制 QQ 图：将观察数据的分位数作为横坐标，理论分位数作为纵坐标，绘制散点图。如果数据符合理论分布，散点图将近似为一条直线。
代码实现
根据上面的理论结果，我们依次进行代码实现。假设数据来自于 gamma 分布。
y <- rgamma(200,10,1)
- 排序数据
sorted_data <- sort(y)
- 计算样本分位数
sample_quantiles <- (1:length(sorted_data) - 0.5) / length(sorted_data)
- 计算理论分位数
首先，对数据进行参数估计，得到估计之后，产生理论分位数。
3.1 参数估计
这里使用极大似然方法估计参数，首先定义对数似然函数：
library(statmod) fn_ig <- function(par= c(mean, dispersion), data) { ll <- -sum(dinvgauss(data, mean= par[1], dispersion= par[2], log= TRUE)) return(ll) }
3.2 参数估计，并计算理论分位数：
pp <- optim(par= c(1,1), fn_ig, data = y) mq <- pp$par[1]; disp <- pp$par[2] theoretical_quantiles <- qinvgauss(sample_quantiles, mean= mq, dispersion= disp)
- 绘制 QQ 图
数据整理并绘制
#数据整理 data_ig = data.frame("x" = theoretical_quantiles, "y" = sorted_data) #绘制 library(ggplot2) ggplot(data_ig,aes(x=x, y=y)) + geom_point(color="#4CA6A3",size=1.2) + geom_abline(intercept = 0, slope = 1,color="#693476",size=1) + # facet_wrap(vars(factor(pc)),scales = "free") + theme_bw() + theme(legend.position = "none") + xlab("Theoretical") + ylab("Sample")
此时，给出了该数据使用逆高斯分布拟合的情况，结果还算可以！读者可以通过一些检验（Kolmogorov–Smirno[3]等），进行进一步检测。
拓展
根据这个思路，小编写了正态和 Gamma 的结果，并通过一个函数将三种分布的整合在一起。最后绘制出三种分布拟合该数据的结果。
还计算了三种分布拟合的 KS 检验的 p 值：
p 值均大于 0.05，接受原假设，即三种分布拟合该数据均可。读者可以将合并的图放到自己论文中，仅需将 y 修改成自己的数据即可！
汇总代码
y <- rgamma(200,10,1) library(statmod) fn_ig <- function(par= c(mean, dispersion), data) { ll <- -sum(dinvgauss(data, mean= par[1], dispersion= par[2], log= TRUE)) return(ll) } fn_normal <- function(par = c(mean, sd), data) { # 正态分布的对数似然函数 ll <- -sum(dnorm(data, mean = par[1], sd = par[2], log = TRUE)) return(ll) } fn_gamma <- function(par = c(shape, rate), data) { # 伽马分布的对数似然函数 ll <- -sum(dgamma(data, shape = par[1], rate = par[2], log = TRUE)) return(ll) } qqdat = function(den = "ig", y = diff_y_all, par= c(1, 1)){ # 排序数据 sorted_data <- sort(y) # 计算样本分位数 sample_quantiles <- (1:length(sorted_data) - 0.5) / length(sorted_data) if(den == "ig"){ # 参数估计 pp <- optim(par= par, fn_ig, data = y) mq <- pp$par[1]; disp <- pp$par[2] # 计算理论分位数 theoretical_quantiles <- qinvgauss(sample_quantiles, mean= mq, dispersion= disp) } else if(den == "norm"){ pp <- optim(par= par, fn_normal, data = y) mq <- pp$par[1]; disp <- pp$par[2] theoretical_quantiles <- qnorm(sample_quantiles, mq, disp) } else if(den == "gamma"){ pp <- optim(par= par, fn_gamma, data = y) mq <- pp$par[1]; disp <- pp$par[2] theoretical_quantiles <- qgamma(sample_quantiles, shape = mq, rate = disp) } data_ig = data.frame("x" = theoretical_quantiles, "y" = sorted_data) return(list(data_ig,"Para"= c(mq,disp))) } # 数据处理 qqplot_den = function(den = "norm",data6 = y,par= c(10, 1)){ u1 = qqdat(den = den, data6, par= par) ig_g = data.frame(u1[[1]]) w1 = ggplot(ig_g,aes(x=x, y=y)) + geom_point(color="#4CA6A3",size=1.2) + geom_abline(intercept = 0, slope = 1,color="#693476",size=1) + theme_bw() + theme(legend.position = "none") + xlab("Theoretical") + ylab("Sample") if(den=="norm"){ kstest = ks.test(data6, "pnorm", mean = u1$Para[1], u1$Para[2])$p.value } else if(den=="ig"){ library(fitdistrplus) kstest = ks.test(data6, "pinvgauss", mean = u1$Para[1], dispersion = u1$Para[2])$p.value } else if(den=="gamma"){ kstest = ks.test(data6, "pgamma", u1$Para[1], u1$Para[2])$p.value } return(list(w1,kstest)) } w1 = qqplot_den(den = "norm",data6 = y,par= c(10, 1)) w2 = qqplot_den(den = "ig",data6 = y,par= c(1, 1)) w3 = qqplot_den(den = "gamma",data6 = y,par= c(1, 1)) cowplot::plot_grid(w1[[1]],w2[[1]],w3[[1]],labels=c("N","IG","G"),nrow=1) # ggsave("qqplot.pdf", width = 8,height = 4)
参考资料
[1]
分位数-分位数图: https://haomenghit.github.io/2019/08/04/Q-Q%E5%9B%BE%E4%BB%8B%E7%BB%8D/
[2]
ggplot2: https://ggplot2.tidyverse.org/
[3]
Kolmogorov–Smirno: https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
