ggplot2中的qqnorm和qqline

Question

假设有一个线性模型LM，我想要画出残差的Q-Q图。通常情况下，我会使用R基础图形：

qqnorm(residuals(LM), ylab="Residuals")
qqline(residuals(LM))

我能够弄清楚如何获得绘图的 qqnorm 部分，但我似乎无法控制 qqline 部分:

ggplot(LM, aes(sample=.resid)) +
    stat_qq()

我怀疑我可能缺少了一些基本的东西，但似乎应该有一种简单的方法可以做到这一点。

编辑：非常感谢下面的解决方案。我稍微修改了代码，从线性模型中提取信息，使得绘图像R基础图形包中的方便绘图一样工作。

ggQQ <- function(LM) # argument: a linear model
{
    y <- quantile(LM$resid[!is.na(LM$resid)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    p <- ggplot(LM, aes(sample=.resid)) +
        stat_qq(alpha = 0.5) +
        geom_abline(slope = slope, intercept = int, color="blue")

    return(p)
}

- Peter

---

注意：对于残差，ggplot需要标准化残差。请参考jlhoward的答案https://dev59.com/oG855IYBdhLWcg3waDiR#19990107。 - qwr

Answer 1

以下代码将会给你想要的图形。ggplot包似乎没有计算qqline参数的代码，所以我不知道是否可以在（可读懂的）一行代码中实现这样的绘图。

qqplot.data <- function (vec) # argument: vector of numbers
{
  # following four lines from base R's qqline()
  y <- quantile(vec[!is.na(vec)], c(0.25, 0.75))
  x <- qnorm(c(0.25, 0.75))
  slope <- diff(y)/diff(x)
  int <- y[1L] - slope * x[1L]

  d <- data.frame(resids = vec)

  ggplot(d, aes(sample = resids)) + stat_qq() + geom_abline(slope = slope, intercept = int)

}

- Aaron

---

完美地运作！我稍微修改了一下代码，直接从线性模型中提取向量。当然，您的解决方案也适用于不是线性模型形式的数据，但我认为其他人可能会需要一个方便的函数来从LM构建qqplot。 - Peter

Answer 2

您还可以使用此功能添加置信区间/置信带（代码的部分内容来自 car:::qqPlot）。

gg_qq <- function(x, distribution = "norm", ..., line.estimate = NULL, conf = 0.95,
                  labels = names(x)){
  q.function <- eval(parse(text = paste0("q", distribution)))
  d.function <- eval(parse(text = paste0("d", distribution)))
  x <- na.omit(x)
  ord <- order(x)
  n <- length(x)
  P <- ppoints(length(x))
  df <- data.frame(ord.x = x[ord], z = q.function(P, ...))

  if(is.null(line.estimate)){
    Q.x <- quantile(df$ord.x, c(0.25, 0.75))
    Q.z <- q.function(c(0.25, 0.75), ...)
    b <- diff(Q.x)/diff(Q.z)
    coef <- c(Q.x[1] - b * Q.z[1], b)
  } else {
    coef <- coef(line.estimate(ord.x ~ z))
  }

  zz <- qnorm(1 - (1 - conf)/2)
  SE <- (coef[2]/d.function(df$z)) * sqrt(P * (1 - P)/n)
  fit.value <- coef[1] + coef[2] * df$z
  df$upper <- fit.value + zz * SE
  df$lower <- fit.value - zz * SE

  if(!is.null(labels)){ 
    df$label <- ifelse(df$ord.x > df$upper | df$ord.x < df$lower, labels[ord],"")
    }

  p <- ggplot(df, aes(x=z, y=ord.x)) +
    geom_point() + 
    geom_abline(intercept = coef[1], slope = coef[2]) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2) 
  if(!is.null(labels)) p <- p + geom_text( aes(label = label))
  print(p)
  coef
}

例子：

Animals2 <- data(Animals2, package = "robustbase")
mod.lm <- lm(log(Animals2$brain) ~ log(Animals2$body))
x <- rstudent(mod.lm)
gg_qq(x)

- Rentrop

---

这非常有帮助。你考虑过在Github上托管你的脚本吗？这样可以很好地引用你的代码。 - Hip Hop Physician

像这样的吗？尽管我不知道为什么你不能引用SO... - Rentrop

非常感谢！我想我说错了一点，我的意思是希望你将它发布在Github上，这样我就可以将其作为R脚本的一部分引入（而不是寻找一种方法来拼接你在Stack Overflow上的帖子）。 - Hip Hop Physician

Answer 3

自3.0版本以来，一个等同于下面的stat_qq_line已经出现在官方ggplot2代码中。

---

旧版本:

从2.0版本开始，ggplot2有一个经过充分记录的扩展接口；因此，我们现在可以轻松地为qqline编写一个新的统计方法（我已经第一次完成了这个任务，所以欢迎改进）：

qq.line <- function(data, qf, na.rm) {
    # from stackoverflow.com/a/4357932/1346276
    q.sample <- quantile(data, c(0.25, 0.75), na.rm = na.rm)
    q.theory <- qf(c(0.25, 0.75))
    slope <- diff(q.sample) / diff(q.theory)
    intercept <- q.sample[1] - slope * q.theory[1]

    list(slope = slope, intercept = intercept)
}

StatQQLine <- ggproto("StatQQLine", Stat,
    # http://docs.ggplot2.org/current/vignettes/extending-ggplot2.html
    # https://github.com/hadley/ggplot2/blob/master/R/stat-qq.r
    
    required_aes = c('sample'),
    
    compute_group = function(data, scales,
                             distribution = stats::qnorm,
                             dparams = list(),
                             na.rm = FALSE) {
        qf <- function(p) do.call(distribution, c(list(p = p), dparams))
        
        n <- length(data$sample)
        theoretical <- qf(stats::ppoints(n))
        qq <- qq.line(data$sample, qf = qf, na.rm = na.rm)
        line <- qq$intercept + theoretical * qq$slope

        data.frame(x = theoretical, y = line)
    } 
)

stat_qqline <- function(mapping = NULL, data = NULL, geom = "line",
                        position = "identity", ...,
                        distribution = stats::qnorm,
                        dparams = list(),
                        na.rm = FALSE,
                        show.legend = NA, 
                        inherit.aes = TRUE) {
    layer(stat = StatQQLine, data = data, mapping = mapping, geom = geom,
          position = position, show.legend = show.legend, inherit.aes = inherit.aes,
          params = list(distribution = distribution,
                        dparams = dparams,
                        na.rm = na.rm, ...))
}

这也可以泛化到分布上（就像stat_qq一样），并可按以下方式使用：

> test.data <- data.frame(sample=rnorm(100, 10, 2)) # normal distribution
> test.data.2 <- data.frame(sample=rt(100, df=2))   # t distribution
> ggplot(test.data, aes(sample=sample)) + stat_qq() + stat_qqline()
> ggplot(test.data.2, aes(sample=sample)) + stat_qq(distribution=qt, dparams=list(df=2)) +
+   stat_qqline(distribution=qt, dparams=list(df=2))

很遗憾，由于qqline在单独的层上，我找不到“重复使用”分布参数的方法，但这只应该是一个小问题。

- phipsgabler

Answer 4

对于线性模型的标准Q-Q诊断，它将标准化残差的分位数绘制为N(0,1)理论量化值。Peter的ggQQ函数绘制了残差。下面的片段对此进行了修改，并添加了一些美化更改，使得绘图更像从plot(lm(...))得到的。

ggQQ = function(lm) {
  # extract standardized residuals from the fit
  d <- data.frame(std.resid = rstandard(lm))
  # calculate 1Q/4Q line
  y <- quantile(d$std.resid[!is.na(d$std.resid)], c(0.25, 0.75))
  x <- qnorm(c(0.25, 0.75))
  slope <- diff(y)/diff(x)
  int <- y[1L] - slope * x[1L]

  p <- ggplot(data=d, aes(sample=std.resid)) +
    stat_qq(shape=1, size=3) +           # open circles
    labs(title="Normal Q-Q",             # plot title
         x="Theoretical Quantiles",      # x-axis label
         y="Standardized Residuals") +   # y-axis label
    geom_abline(slope = slope, intercept = int, linetype="dashed")  # dashed reference line
  return(p)
}

使用示例：

# sample data (y = x + N(0,1), x in [1,100])
df <- data.frame(cbind(x=c(1:100),y=c(1:100+rnorm(100))))
ggQQ(lm(y~x,data=df))

- jlhoward

Answer 5

使用最新的ggplot2版本（>=3.0），实现了新函数stat_qq_line (https://github.com/tidyverse/ggplot2/blob/master/NEWS.md)，可以通过以下方式添加模型残差的qq线：

library(ggplot2)
model <- lm(mpg ~ wt, data=mtcars)
ggplot(model, aes(sample = rstandard(model))) + geom_qq() + stat_qq_line()

rstandard(model)用于获取标准化残差。(感谢@jlhoward和@qwr)

如果出现'Error in stat_qq_line() : could not find function "stat_qq_line"'的错误提示，说明您的ggplot2版本太旧了，可以通过升级ggplot2包来解决：install.packages("ggplot2")。

- LmW.

---

现在ggplot2 3.0.0版本已经稳定发布，因此您现在可以从CRAN版本中获取此功能。 - Droplet

Ggplot期望如jlhoward所描述的标准化残差。因此，请使用rstandard(model)而不是.resid。 - qwr

Answer 6

为什么不使用以下内容？

假设有一个向量，比如说，

myresiduals <- rnorm(100) ^ 2

ggplot(data=as.data.frame(qqnorm( myresiduals , plot=F)), mapping=aes(x=x, y=y)) + 
    geom_point() + geom_smooth(method="lm", se=FALSE)

但是使用传统的图形函数来支持ggplot2似乎很奇怪。

我们不能通过从我们想要量化图的向量开始，然后在ggplot2中应用相应的“stat”和“geom”函数来获得同样的效果吗？

Hadley Wickham是否监控这些帖子？也许他可以向我们展示更好的方法。

- Jacob Wegelin

---

散点图类似于qqnorm()的Q-Q图，但geom_smooth添加的线与qqline()给出的线不同。另一方面，Aaron和@jlhoward提供的解决方案给出了类似于基本R的图形。您能否评论一下是否是我的数据导致它表现不佳？ - ktyagi

Answer 7

您可以学习老一辈使用普通概率纸的方法。仔细观察 ggplot()+stat_qq() 图形，可以使用 geom_abline() 添加参考线，就像这样：

df <- data.frame( y=rpois(100, 4) )

ggplot(df, aes(sample=y)) +
  stat_qq() +
  geom_abline(intercept=mean(df$y), slope = sd(df$y))

- Mike Anderson

---

一个样本与理论量化值的Q-Q图不应该有参考线y=x吗？ - qwr

Answer 8

ggplot2 v.3.0.0现在具有qqline stat。来自帮助页面：

df <- data.frame(y = rt(200, df = 5))
p <- ggplot(df, aes(sample = y))
p + stat_qq() + stat_qq_line()

!ggplot2 v3.0.0示例，相当于qqnorm加上abline的统计图示例]1

- Richard Careaga