This function plots a QQplot as ggplot in R - Parts of the code copied from `car:::qqPlot`

gg_qq.r

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
}

test.R

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

# See http://stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2/#27191036

I'd like to include this in my package (userfriendlyscience, it has a 'normalityAssessment' function which is a wrapper for useful functions to, well, assess normality :-)). Would that be ok, and if so, how would you prefer to be credited?

Sure! You are free to use it. Just cite my github-account. That is enough.

Great, thank you very much! Also on behalf of the researchers who will hopefully benefit from your function in the future! 😄

You might want to add the ellipsis to the call of the d-function, too.

line 21:
SE <- (coef[2]/d.function(df$z)) * sqrt(P * (1 - P)/n)

to
SE <- (coef[2]/d.function(df$z, ...)) * sqrt(P * (1 - P)/n)

Currently, either the evaluation fails because of missing arguments or the confidence intervals are messed up for distributions like weibull, gamma, or lnorm.

These confidence intervals are based on the approach of John Fox in Applied Regression Analysis book?