Сложный регрессионный график в R

10

Мне нужно нарисовать сложную графику для визуального анализа данных. У меня есть 2 переменные и большое количество случаев (> 1000). Например (число равно 100, если дисперсия меньше "нормальной"):

x <- rnorm(100,mean=95,sd=50)
y <- rnorm(100,mean=35,sd=20)
d <- data.frame(x=x,y=y)

1) Мне нужно построить исходные данные с размером точки, соответствующей относительной частоте совпадений, поэтому plot(x,y)это не вариант - мне нужны размеры точек. Что нужно сделать для этого?

2) На том же графике мне нужно построить эллипс с доверительным интервалом 95% и линию, представляющую изменение корреляции (не знаю, как правильно назвать ее) - что-то вроде этого:

library(corrgram)
corrgram(d, order=TRUE, lower.panel=panel.ellipse, upper.panel=panel.pts)

correlogramm

но с обоими графиками на одном участке.

3) Наконец, мне нужно нарисовать результирующую модель линейной регрессии поверх всего этого:

r<-lm(y~x, data=d)
abline(r,col=2,lwd=2)

но с диапазоном ошибок ... как на QQ-plot:

QQ-график

но для ошибок подгонки, если это возможно.

Итак, вопрос:

Как добиться всего этого на одном графике?

Юрий Петровский
источник

Ответы:

29

Похоже ли изображение ниже на то, чего вы хотите достичь?

введите описание изображения здесь

Вот обновленный код R, после ваших комментариев:

do.it <- function(df, type="confidence", ...) {
  require(ellipse)
  lm0 <- lm(y ~ x, data=df)
  xc <- with(df, xyTable(x, y))
  df.new <- data.frame(x=seq(min(df$x), max(df$x), 0.1))
  pred.ulb <- predict(lm0, df.new, interval=type)
  pred.lo <- predict(loess(y ~ x, data=df), df.new)
  plot(xc$x, xc$y, cex=xc$number*2/3, xlab="x", ylab="y", ...)
  abline(lm0, col="red")
  lines(df.new$x, pred.lo, col="green", lwd=1.5)
  lines(df.new$x, pred.ulb[,"lwr"], lty=2, col="red")
  lines(df.new$x, pred.ulb[,"upr"], lty=2, col="red")    
  lines(ellipse(cor(df$x, df$y), scale=c(sd(df$x),sd(df$y)), 
        centre=c(mean(df$x),mean(df$y))), lwd=1.5, col="green")
  invisible(lm0)
}

set.seed(101)
n <- 1000
x <- rnorm(n, mean=2)
y <- 1.5 + 0.4*x + rnorm(n)
df <- data.frame(x=x, y=y)

# take a bootstrap sample
df <- df[sample(nrow(df), nrow(df), rep=TRUE),]

do.it(df, pch=19, col=rgb(0,0,.7,.5))

И вот версия ggplotized

введите описание изображения здесь

производится со следующим фрагментом кода:

xc <- with(df, xyTable(x, y))
df2 <- cbind.data.frame(x=xc$x, y=xc$y, n=xc$number)
df.ell <- as.data.frame(with(df, ellipse(cor(x, y), 
                                         scale=c(sd(x),sd(y)), 
                                         centre=c(mean(x),mean(y)))))
library(ggplot2)

ggplot(data=df2, aes(x=x, y=y)) + 
  geom_point(aes(size=n), alpha=.6) + 
  stat_smooth(data=df, method="loess", se=FALSE, color="green") + 
  stat_smooth(data=df, method="lm") +
  geom_path(data=df.ell, colour="green", size=1.2)

Его можно настроить немного больше, добавив индексы соответствия модели, такие как расстояние Кука, с эффектом затенения цветом.

хл
источник
1
@chl +1, хороший график и короткий код.
mpiktas
@mpiktas Спасибо. Это привело меня к пониманию, что я не работал с правильным образцом, на самом деле :-)
chl
Почти выглядит как тот, который мне нужен, но с реальными числами я столкнулся со следующими проблемами: 1) df.new <- data.frame(x = seq(min(x), max(x), 0.1))лучше. 2) Эллипс рисуется в позиции 0; 0, что неверно, и это s size is also strange (too small). Also tryed библиотека (автомобиль) dataEllipse (df y, уровни = 0,95: 1, lty = 2) `, но он отбрасывает все. 3) Кривая (как на корлограмме) отсутствует. Я почти воспроизвел это по телефону, но диапазон данных неверен. Используйте первые 2 строки из моего кода вместо вашей для воспроизведения. Икс,dеlibrary(car) cr.plots(m0)
Юрий Петровский
@Yuriy Хорошо, я обновлю свой код (пока что не нужно вносить какие-либо изменения), но я не могу понять, как мы можем перекрываться с действительными значениями случайных величин с вашей настройкой ; это причина, почему я использую Boostrap с заменой (это гарантирует, что ~ 2/3 оригинальных единиц присутствуют). действительно предоставляет те же возможности, что и в пакете, но, вероятно, это не так легко настроить. Я полагаю, что наложенная кривая - это всего лишь лесс , поэтому добавить ее нетрудно. (Икс,Y)car::dataEllipseellipse
chl
2
@Tal Интерпретация эллипса такая же, как и в corrgramпакете: она показывает 95% парную доверительную область, предполагая двумерное нормальное распределение с центром в среднем и масштабированное с помощью SD (x) и SD (y). Я не большой поклонник этого, когда используется на графике рассеяния, хотя. Но см. Murdoch & Chow, Графическое отображение больших корреляционных матриц , Am Stat (1996) 50: 178 или Friendly, Corrgrams: Исследовательские дисплеи для корреляционных матриц , Am Stat (2002) 56: 316.
ЧЛ
2

Для точки 1 просто используйте cexпараметр на графике, чтобы установить размер точки.

Например

x = rnorm(100)
plot(x, pch=20, cex=abs(x))

Чтобы иметь несколько графиков на одном графике, используйте par(mfrow=c(numrows, numcols))равномерно распределенный макет или layoutсделайте более сложные.

Nico
источник
1
+1 за подсказку cex, но я думаю, что ОП хочет все вещи в одном и том же регионе, а не в отдельных.
ЧЛ
Ааа ... теперь я понимаю вопрос. Ну, тогда он может просто использовать curveили pointsпересекать три графика;)
Нико