在ggplot2中的qqnorm和qqline
说有一个线性模型LM,我想要一个残差的qq图。 通常我会使用R基本graphics:
qqnorm(residuals(LM), ylab="Residuals") qqline(residuals(LM))
我可以弄清楚如何获得情节的qqnorm部分,但我似乎无法pipe理qqline:
ggplot(LM, aes(sample=.resid)) + stat_qq()
我怀疑我错过了一些非常基本的东西,但似乎应该有一个简单的方法来做到这一点。
编辑:非常感谢下面的解决scheme。 我已经修改了代码(非常轻微地)从线性模型中提取信息,以使绘图像R基本graphics包中的便捷图一样工作。
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) }
下面的代码会给你你想要的阴谋。 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) }
您也可以使用此function添加置信区间/置信区间(从car:::qqPlot
复制的代码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)){ Qx <- quantile(df$ord.x, c(0.25, 0.75)) Qz <- q.function(c(0.25, 0.75), ...) b <- diff(Qx)/diff(Qz) coef <- c(Qx[1] - b * Qz[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)
线性模型的标准QQ诊断标绘出标准化残差与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 QQ", # 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))
为什么不是以下?
给定一些vector,比如说,
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)
但似乎很奇怪,我们必须使用传统的graphicsfunction来支撑ggplot2。
我们不能通过从我们想要分位点图的vector开始,然后在ggplot2中应用相应的“统计”和“几何”函数来得到相同的效果吗?
Hadley Wickham是否监视这些post? 也许他可以给我们一个更好的方法。
从版本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是在一个单独的层次上,我找不到一个“重用”分布参数的方法,但这应该只是一个小问题)。