Domanda

Di 'hanno un modello LM lineare che voglio una trama qq dei residui. Normalmente avrei usare la grafica di base R:

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

riesco a capire come ottenere la parte qqnorm della trama, ma io non riesco a gestire il qqline:

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

ho il sospetto che mi manca qualcosa piuttosto semplice, ma sembra che ci dovrebbe essere un modo semplice di fare questo.

EDIT: Molte grazie per la soluzione qui di seguito. Ho modificato il codice (leggermente) per estrarre le informazioni dal modello lineare in modo che la trama funziona come la trama convenienza nel pacchetto grafico di base 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)
}
È stato utile?

Soluzione

Il seguente codice vi darà la trama che si desidera. Il pacchetto ggplot non sembra contenere codice per il calcolo dei parametri del qqline, quindi non so se è possibile ottenere un tale trama in un (comprensibile) one-liner.

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)

}

Altri suggerimenti

È inoltre possibile aggiungere gli intervalli di confidenza / bande di confidenza con questa funzione (Parti del codice copiato da 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
}

Esempio:

Animals2 <- data(Animals2, package = "robustbase")
mod.lm <- lm(log(Animals2$brain) ~ log(Animals2$body))
x <- rstudent(mod.lm)
gg_qq(x)

entrare descrizione dell'immagine qui

Il Q-Q standard diagnostico per lineari modelli appezzamenti quantili delle standardizzato residui vs quantili teorici di N (0,1). @ Funzione ggQQ di Peter Traccia i residui. Il frammento di seguito ammenda che e aggiunge alcuni cosmetici cambia per rendere la trama più simile a quello che si potrebbe ottenere da 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)
}

Esempio di utilizzo:

# 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))

A partire dalla versione 2.0, ggplot2 ha un'interfaccia ben documentata per l'estensione; così ora possiamo facilmente scrivere un nuovo stat per il qqline da solo (che ho fatto per la prima volta, quindi i miglioramenti sono benvenuto ):

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, ...))
}

Questo generalizza anche sulla distribuzione (esattamente come stat_qq fa), e può essere utilizzato come segue:

> 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))

(Purtroppo, dal momento che il qqline è su un livello separato, non ho potuto trovare un modo per "riuso" i parametri di distribuzione, ma che dovrebbe essere solo un problema minore.)

Perché non il seguente?

Dato un vettore, per esempio,

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)

Ma sembra strano che dobbiamo utilizzare una grafica tradizionale funzionare per puntellare ggplot2.

Non possiamo ottenere lo stesso effetto in qualche modo iniziando con il vettore per il quale vogliamo che la trama quantile e quindi applicando il "stat" appropriata e funzioni "geom" in ggplot2?

non Hadley Wickham monitorare questi posti? Forse ci può mostrare un modo migliore.

Con l'ultima versione ggplot2 (> = 3.0), nuova funzione stat_qq_line è implementata ( https://github.com/tidyverse/ggplot2/blob/master/NEWS.md ) e una linea QQ per residui modello possono essere aggiunti con:

library(ggplot2)
model <- lm(mpg ~ wt, data=mtcars)
ggplot(model, aes(sample = rstandard(model))) + geom_qq() + stat_qq_line()
è necessario

rstandard(model) per ottenere il residuo standardizzato. (@Jlhoward credito e @qwr)

Se si ottiene un 'errore in stat_qq_line (): non riusciva a trovare la funzione 'stat_qq_line'', la versione ggplot2 è troppo vecchio e si può risolvere il problema aggiornando il proprio pacchetto di ggplot2:. install.packages("ggplot2")

Si potrebbe rubare una pagina dai veterani che hanno fatto questa roba con la carta di probabilità normale. Uno sguardo attento alla ggplot () + stat_qq () grafico suggerisce che una linea di riferimento può essere aggiunta con geom_abline (), come questo

df <- data.frame( y=rpois(100, 4) )

ggplot(df, aes(sample=y)) +
  stat_qq() +
  geom_abline(intercept=mean(df$y), slope = sd(df$y))

v.3.0.0 ggplot2 ora ha una statistica qqline. Dalla pagina di aiuto:

df <- data.frame(y = rt(200, df = 5))
p <- ggplot(df, aes(sample = y))
p + stat_qq() + stat_qq_line()

! Ggplot2 v3.0.0 Esempio statistiche equivalente a qqnorm più abline ] 1

Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top