ggplot2のqqnormとqqline
https://stackoverflow.com/questions/4357031
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
残差のQQプロットが必要な線形モデルlmを持っているとします。通常、Rベースグラフィックスを使用します。
qqnorm(residuals(LM), ylab="Residuals")
qqline(residuals(LM))
プロットのqqnorm部分を取得する方法を理解できますが、qqlineを管理することはできません。
ggplot(LM, aes(sample=.resid)) +
stat_qq()
私はかなり基本的なものが欠けているのではないかと思いますが、これを行う簡単な方法があるべきだと思われます。
編集： 以下の解決策に感謝します。コードを（非常にわずかに）変更して、線形モデルから情報を抽出して、プロットがRベースグラフィックスパッケージのコンビニエンスプロットのように機能するようにしました。
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)
}
他のヒント
また、この関数で信頼区間/信頼帯を追加することもできます（コードの一部がコピーされました 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)
線形モデルの標準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 Q-Q", # 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))
バージョン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は別のレイヤーにあるため、分布パラメーターを「再利用」する方法を見つけることができませんでしたが、それはわずかな問題にすぎません。）
なぜ次のことではないのですか？
いくつかのベクトルを与えられた、たとえば、
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)
しかし、従来のグラフィック機能を使用してGGPLOT2を支えなければならないのは奇妙に思えます。
ギャンシタルプロットが必要なベクトルから始めて、GGPLOT2に適切な「統計」関数と「GEOM」関数を適用することで、どういうわけか同じ効果を得ることができませんか？
ハドリー・ウィッカムはこれらの投稿を監視していますか？たぶん彼は私たちにもっと良い方法を示すことができます。
最新のggplot2バージョン（> = 3.0）で、新しい機能 stat_qq_line 実装されています（https://github.com/tidyverse/ggplot2/blob/master/news.md）およびモデル残差のQQラインを次のように追加できます。
library(ggplot2)
model <- lm(mpg ~ wt, data=mtcars)
ggplot(model, aes(sample = rstandard(model))) + geom_qq() + stat_qq_line()
rstandard(model) 標準化された残差を取得するには必要です。 （クレジット@jlhowardと@qwr）
stat_qq_line（）で 'エラーが発生した場合：function "stat_qq_line"'を見つけることができなかった場合、ggplot2バージョンが古すぎて、ggplot2パッケージをアップグレードして修正できます。 install.packages("ggplot2") .
通常の確率論文でこのようなことをした昔の人からページを盗むことができます。 ggplot（）+stat_qq（）グラフィックを慎重に見ると、このようなgeom_abline（）で参照線を追加できることが示唆されています。
df <- data.frame( y=rpois(100, 4) )
ggplot(df, aes(sample=y)) +
stat_qq() +
geom_abline(intercept=mean(df$y), slope = sd(df$y))
ggplot2 v.3.0.0にはqqline統計があります。ヘルプページから：
df <- data.frame(y = rt(200, df = 5))
p <- ggplot(df, aes(sample = y))
p + stat_qq() + stat_qq_line()
