Skip to content

Instantly share code, notes, and snippets.

@rentrop
Last active December 3, 2023 06:01
  • Star 8 You must be signed in to star a gist
  • Fork 3 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save rentrop/d39a8406ad8af2a1066c to your computer and use it in GitHub Desktop.
This function plots a QQplot as ggplot in R - Parts of the code copied from `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)
# See http://stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2/#27191036
@rentrop
Copy link
Author

rentrop commented Sep 26, 2015

Parts of the code copied from car:::qqPlot
This is a way to plot qqnorm and qqline (including confidence intervals) using ggplot2.

For other ways to do it see http://stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2/

@Matherion
Copy link

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?

@rentrop
Copy link
Author

rentrop commented Mar 19, 2016

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

@Matherion
Copy link

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

@lkreft
Copy link

lkreft commented Mar 20, 2017

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.

@bravo-diego
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment