2010-12-05 5 views
48

残差のqqプロットが必要な線形モデルLMがあります。通常、私はRベースのグラフィックスを使用します。qqnormとqqline in ggplot2

qqnorm(residuals(LM), ylab="Residuals") 
qqline(residuals(LM)) 

私はプロットのqqnormの一部を取得する方法を見つけ出すことができますが、私はqqlineを管理するように見えることはできません。

ggplot(LM, aes(sample=.resid)) + 
    stat_qq() 

私は疑います私はかなり基本的な何かを欠いているが、これを行う簡単な方法であるべきであるように思える。

EDIT:下記のソリューションに感謝します。私は線形モデルから情報を抽出するようにコードを修正しました(非常にわずかに)ので、プロットは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) 
} 

答えて

50

次のコードは、必要なプロットを提供します。 ggplotパッケージはqqlineのパラメータを計算するためのコードを含んでいないようですので、(理解しやすい)1ライナーでそのようなプロットを達成することが可能かどうかはわかりません。

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) 

} 
+0

完璧に動作します!私は線形モデルからベクトルを直接抽出するコードをわずかに修正する自由を取った。もちろん、あなたのソリューションは線形モデルの形ではないデータでも動作しますが、他の誰かがLMからqqplotを構築するための便利な機能を望むと思っています。 – Peter

9

以下の理由はありますか?いくつかのベクトルを考えると

、言う、

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で適切な "stat"関数と "geom"関数を適用すると、同じ効果が得られませんか?

ハドリーウィッカムはこれらの投稿を監視していますか?たぶん彼は私たちにもっと良い方法を教えてくれるかもしれません。

+0

散布図はqqnorm()のq-qプロットに似ていますが、geom_smoothによって追加された行はqqline()によって与えられた行と同じではありません。一方、Aaronと@jlhowardによって与えられた解は、基底のRに類似したプロットを与える。それが自分のデータであるとコメントすることができます。なぜなら、それは誤動作しているからです。 – ktyagi

10

線形モデルの標準Q-Q診断は、の標準化された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)) 
19

また、この機能で信頼区間/信頼帯を追加することができます(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) 

enter image description here

+1

これは非常に便利です。 Githubであなたのスクリプトをホストすることについて考えましたか?コードを適切に配置するといいでしょう。 –

+1

https://gist.github.com/rentrop/d39a8406ad8af2a1066c?私はなぜあなたがサイトを見つけることができないのか知りません。 – Rentrop

+3

ありがとう!私はちょっとミスをしたと思うが、Githubに投稿されているのはいいと思っていたので、Stack Overflowの投稿をスプライスする方法を見つけるのではなく、Rスクリプトの一部として取り込むことができる。 –

11

バージョン2以降です。0の場合、ggplot2には拡張のための十分に文書化されたインタフェースがあります。私たちは今、簡単に(改善がwelcomeているので、私が最初にやったもの)自体によってqqlineのための新たなSTATを書き込むことができます。

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は別のレイヤー上にあるため、残念ながら、私は分布パラメータを「再利用」する方法を見つけることができませんでしたが、それが唯一のマイナーな問題でなければなりません。)

1

あなたはページを盗むことができます普通の確率紙でこれをやった昔の人。 gogplot()+ stat_qq()グラフィックを注意深く見ると、このようにgeom_abline()で参照線を追加できることがわかります。