R语言画直线加置信区间 r语言 置信区间图

前言

检验数据的分布有很多种方法，如卡方检验、K-S 检验等。

而以图形的角度来说，可以使用 Q-Q 图或 P-P 图来检验数据是否服从：

• beta 分布：
• t 分布
• 卡方分布
• gamma 分布
• 正态分布
• 均匀分布
• logistic 分布等

Q-Q 图通过将两个概率分布的相同分位数点的值映射为 x 和 y 轴。

如果两个分布比较相似，那么图上的点大致落在 y=x 直线上；如果这两个分布线性相关，则点大致分布在一条直线上，但不一定在 y=x 上。

如果样本数据近似服从正态分布，则图上的点大致落在一条直线上，该直线的斜率为标准差，截距为均值。

P-P 图根据变量的累积比例与特定分布的累积比例之间的关系来绘制图形，可以检验数据是否服从指定的分布。

如果数据服从指定分布，则图中的点近似地呈一条直线，如果不呈直线，但是存在某种规律，那么可以通过某种变换近似的逼近指定分布

Q-Q 图和 P-P 图除了检验方法不同之外，其他完全一样

示例

Q-Q 图

1. ggplot2

可以使用 ggplot2 的 geom_qq 和 geom_qq_line 绘制 Q-Q 图

tibble(y = rt(200, df = 5)) %>%
  ggplot(aes(sample = y)) +
  geom_qq() + geom_qq_line()

使用 MASS 包的 fitdistr 函数来估计分布的参数，并传递到理论分布中

df <- tibble(y = rt(200, df = 5))

params <- as.list(MASS::fitdistr(df$y, "t")$estimate)

ggplot(df, aes(sample = y)) +
  geom_qq(distribution = qt, dparams = params["df"]) +
  geom_qq_line(distribution = qt, dparams = params["df"])

我们可以看到，直线与点之间的趋势更加紧密

我们换一个分布，就 beta 分布吧

df <- tibble(y = rbeta(200, 1, 2))

ggplot(df, aes(y)) +
  geom_density()

params <- as.list(MASS::fitdistr(df$y, "beta", 
                                 start=list(shape1=1, shape2=2))$estimate)

ggplot(df, aes(sample = y)) +
  geom_qq(distribution = qbeta, dparams = params) +
  geom_qq_line(distribution = qbeta, dparams = params)

绘制分组数据的 Q-Q 图

ggplot(mtcars, aes(sample = mpg, colour = factor(cyl))) +
  stat_qq() +
  stat_qq_line()

2. qqplotr

我们也可以使用扩展包 qqplotr 来绘制 Q-Q 图

if (!require(qqplotr)) {
  install.packages("qqplotr")
}
require(qqplotr)

df <- tibble(x = rnorm(200))

ggplot(df, aes(sample = x)) +
  stat_qq_band() +
  stat_qq_line() +
  stat_qq_point() +
  labs(x = "Theoretical", y = "Sample")

其中 stat_qq_band 用于绘制置信区间，有如下几种方法:

• "pointwise"：使用正态分布构造
• "boot"：基于 boostrap 构造
• "ks"：基于 Kolmogorov-Smirnov 检验
• "ts"：构造 tail-sensitive 的置信区间

我们可以使用 bandType 参数来指定

ggplot(df, aes(sample = x)) +
  geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.7) +
  geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.7) +
  geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.7) +
  geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.7) +
  stat_qq_line() +
  stat_qq_point() +
  scale_fill_manual(values = c("#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3")) +
  labs(x = "Theoretical", y = "Sample")

我们也可以使用其他理论分布构建 Q-Q 图，分布的参数作为一个 list 传递给 dparams 参数，这一点与 ggplot2 是一样的

例如，指数分布

df <- tibble(x = rexp(300, 2))

dist <- "exp"
params <- list(rate = 2)

ggplot(df, aes(sample = x)) +
  stat_qq_band(distribution = dist, dparams = params, fill = "#8da0cb") +
  stat_qq_line(distribution = dist, dparams = params) +
  stat_qq_point(distribution = dist, dparams = params) +
  labs(x = "Theoretical", y = "Sample")

qqplotr 还可以将参考线水平化，以减少视觉上的偏差，可以设置 detrend = TRUE

ggplot(df, aes(sample = x)) +
  stat_qq_band(distribution = dist, dparams = params, fill = "#8da0cb", detrend = TRUE) +
  stat_qq_line(distribution = dist, dparams = params, detrend = TRUE) +
  stat_qq_point(distribution = dist, dparams = params, detrend = TRUE) +
  labs(x = "Theoretical", y = "Sample")

还支持分面操作

dist <- "beta"
params <- list(shape1 = 1, shape2 = 2)

ggplot(diamonds, aes(sample = price, colour = cut, fill = cut)) +
  stat_qq_band(distribution = dist, dparams = params, alpha = 0.6) +
  stat_qq_point(distribution = dist, dparams = params) +
  stat_qq_line(distribution = dist, dparams = params) +
  facet_wrap(vars(cut)) +
  labs(x = "Theoretical", y = "Sample")

P-P 图

我们还是使用 qqplotr 来绘制 P-P 图

例如，检验数据是否服从标准正态分布

df <- tibble(x = rnorm(200))

ggplot(df, aes(sample = x)) +
  stat_pp_band() +
  stat_pp_line() +
  stat_pp_point() +
  labs(x = "Probability Points", y = "Cumulative Probability")

df <- tibble(x = rnorm(200, mean = 1, sd = 2))

params <- list(mean = 5, sd = 2)

ggplot(df, aes(sample = x)) +
  stat_pp_band(dparams = params) +
  stat_pp_line() +
  stat_pp_point(dparams = params) +
  labs(x = "Probability Points", y = "Cumulative Probability")

注意：stat_pp_line 不支持传入分布，但是可以设置 ab 参数的值来指定直线的斜率和截距

ggplot(df, aes(sample = x)) +
  stat_pp_band(dparams = params) +
  stat_pp_line(ab = c(0, 0.1)) + # c(截距,斜率)
  stat_pp_point(dparams = params) +
  labs(x = "Probability Points", y = "Cumulative Probability")

检验其他分布

df <- tibble(x = rgamma(300, shape = 1, rate = 0.2))

dist <- "gamma"
params <- list(shape = 1, rate = 0.2)

ggplot(df, aes(sample = x)) +
  stat_pp_band(distribution = dist, dparams = params, fill = "#8da0cb", detrend = TRUE) +
  stat_pp_line(detrend = TRUE, colour = "#66c2a5") +
  stat_pp_point(distribution = dist, dparams = params, colour = "#e78ac3", detrend = TRUE) +
  ylim(c(-0.5, 0.5)) +
  labs(x = "Probability Points", y = "Cumulative Probability")