-
-
Save rentrop/d39a8406ad8af2a1066c to your computer and use it in GitHub Desktop.
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 |
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?
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/