qqnorm et qqline dans ggplot2
Question
Say ont un modèle linéaire LM que je veux un terrain qq des résidus. Normalement, j'utiliser les graphiques de base R:
qqnorm(residuals(LM), ylab="Residuals")
qqline(residuals(LM))
Je peux comprendre comment obtenir la partie qqnorm de la parcelle, mais je ne peux pas sembler gérer le qqline:
ggplot(LM, aes(sample=.resid)) +
stat_qq()
Je pense que je me manque quelque chose assez basique, mais il semble que il devrait y avoir un moyen facile de le faire.
EDIT: Merci beaucoup pour la solution ci-dessous. J'ai modifié le code (très légèrement) pour extraire les informations du modèle linéaire afin que l'intrigue fonctionne comme le tracé de la commodité dans le package graphique de 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)
}
La solution
Le code suivant vous donnera le tracé que vous voulez. Le paquet ggplot ne semble pas contenir du code pour le calcul des paramètres du qqline, donc je ne sais pas s'il est possible de réaliser un tel complot dans un (compréhensible) en une ligne.
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)
}
Autres conseils
Vous pouvez également ajouter Intervalles de confiance / bandes de confiance avec cette fonction (parties du code copié de 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
}
Exemple:
Animals2 <- data(Animals2, package = "robustbase")
mod.lm <- lm(log(Animals2$brain) ~ log(Animals2$body))
x <- rstudent(mod.lm)
gg_qq(x)
Le diagnostic standard Q-Q pour les modèles linéaires terrains quantiles des résidus normalisé par rapport aux quantiles théoriques de N (0,1). @ La fonction de Peter ggQQ trace les résidus. L'extrait ci-dessous amende honorable et ajoute que quelques changements cosmétiques pour rendre l'intrigue plus comme ce que l'on pourrait obtenir de 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)
}
Exemple d'utilisation:
# 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))
Depuis la version 2.0, ggplot2 dispose d'une interface bien documentée pour l'extension; de sorte que nous pouvons maintenant facilement écrire une nouvelle stat pour la qqline par lui-même (que je l'ai fait pour la première fois, des améliorations sont ):
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, ...))
}
généralise sur la distribution (exactement comme stat_qq
fait), et peut être utilisé comme suit:
> 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))
(Malheureusement, depuis le qqline est sur une couche séparée, je ne pouvais pas trouver un moyen de « réutilisation » les paramètres de distribution, mais qui ne devrait être un problème mineur.)
Pourquoi pas ce qui suit?
Étant donné un vecteur, disons,
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)
Mais il semble étrange que nous devons utiliser une fonction graphiques traditionnels pour soutenir ggplot2.
On ne peut pas obtenir le même effet en quelque sorte en commençant par le vecteur pour lequel nous voulons que le diagramme quantile puis appliquer le appropriée « stat » et les fonctions « geom » dans ggplot2?
Est-ce que Hadley Wickham surveiller ces messages? Peut-être qu'il peut nous montrer une façon meilleure.
Avec la dernière version ggplot2 (> = 3.0), nouvelle Vous pourriez voler une page des anciens qui ont fait ce genre de choses avec du papier de probabilité normale. Un regard attentif à un ggplot () + stat_qq () graphique indique qu'une ligne de référence peut être ajouté avec geom_abline (), comme celui-ci stat_qq_line
de fonction est implémentée (
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 a maintenant une stat qqline. De la page d'aide:
df <- data.frame(y = rt(200, df = 5))
p <- ggplot(df, aes(sample = y))
p + stat_qq() + stat_qq_line()
! Ggplot2 v3.0.0 Exemple stats équivalent à qqnorm, plus abline ] 1