统计问题｜绘制任意分布的 QQ 图
点击下方公众号，回复资料分享，收获惊喜
简介
最近科研论文的审稿意见，需要对数据的拟合情况进行说明。分位数-分位数图[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)
参考资料
分位数-分位数图: 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
相关推荐
colourpicker包：图形颜色拾取器
visdat包：助你一眼看穿数据结构和缺失值！
report包：助你自动出统计报告！
reticulate包｜数据科学者的福音