ggplot2 quantile-quantile plot that contains confidence bounds

geom_qq_conf.R

# code copied and amended from http://stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2/

library(ggplot2)

gg_qq_conf <- 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_line(aes(y = lower), lty=2) +
    geom_line(aes(y = upper), lty=2)+
  #  geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2) +
    xlab("Theoretical quantiles") +
    ylab("Sample quantiles")
  if(!is.null(labels)) p <- p + geom_text( aes(label = label))
  return(p)
  #coef
}

data(Animals)
mod.lm <- lm(data=Animals,
             log(brain) ~ log(body))

x <- rstudent(mod.lm) # studentised residuals

gg_qq_conf(x, labels=NULL) + theme_bw() + coord_equal()