Wykres regresji złożonej w R.

10

Muszę narysować złożoną grafikę do wizualnej analizy danych. Mam 2 zmienne i dużą liczbę przypadków (> 1000). Na przykład (liczba wynosi 100, jeśli dyspersja jest mniej „normalna”):

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

1) Muszę wykreślić surowe dane z rozmiarem punktu, odpowiadającym względnej częstotliwości zbieżności, więc plot(x,y)nie ma opcji - potrzebuję rozmiarów punktów. Co należy zrobić, aby to osiągnąć?

2) Na tym samym wykresie muszę wykreślić elipsę 95% przedziału ufności i linię reprezentującą zmianę korelacji (nie wiem jak poprawnie ją nazwać) - coś takiego:

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

correlogramm

ale z dwoma wykresami na jednym wykresie.

3) Na koniec muszę narysować wynikowy model regresji liniowej nad tym wszystkim:

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

ale z zakresem błędów ... coś jak na wykresie QQ:

Wykres QQ

ale w przypadku błędów dopasowania, jeśli jest to możliwe.

Pytanie brzmi:

Jak to wszystko osiągnąć na jednym wykresie?

Jurij Pietrowski
źródło

Odpowiedzi:

29

Czy obrazek poniżej wygląda tak, jak chcesz?

wprowadź opis zdjęcia tutaj

Oto zaktualizowany kod R z następującymi komentarzami:

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

A oto wersja ggplotized

wprowadź opis zdjęcia tutaj

wyprodukowany z następującym fragmentem kodu:

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)

Można go jeszcze bardziej dostosować, dodając wskaźniki dopasowania modelu, takie jak odległość Cooka, z efektem cieniowania kolorów.

chl
źródło
1
@chl +1, ładny wykres i krótki kod.
mpiktas,
@mpiktas Thanks. Doprowadziło mnie to do wniosku, że nie pracowałem z odpowiednią próbką, w rzeczywistości :-)
chl
df.new <- data.frame(x = seq(min(x), max(x), 0.1))s size is also strange (too small). Also tryed x,refalibrary(car) cr.plots(m0)
(x,y)car::dataEllipseellipse
2
@Tal Interpretacja elipsy jest taka sama jak w corrgrampakiecie: pokazuje ona 95% obszar ufności parami, zakładając dwuwymiarowy rozkład normalny wyśrodkowany na średniej i skalowany przez SD (x) i SD (y). Jednak nie jestem wielkim fanem tego, gdy używa się go na wykresie rozrzutu. Ale patrz Murdoch i Chow, Graficzny wyświetlacz dużych macierzy korelacji , Am Stat (1996) 50: 178 lub Friendly, Corrgrams: Wyświetlacze eksploracyjne dla macierzy korelacji , Am Stat (2002) 56: 316.
chl
2

Dla punktu 1 wystarczy użyć cexparametru na wykresie, aby ustawić rozmiar punktu.

Na przykład

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

Aby mieć wiele wykresów na jednym wykresie, użyj par(mfrow=c(numrows, numcols))równomiernie rozmieszczonego układu lub layoutwykonaj bardziej złożone.

Nico
źródło
1
+1 za wskazówkę cex, ale myślę, że OP chce wszystkich rzeczy na tym samym obszarze kreślenia, a nie na osobnych.
chl
Ahh ... teraz rozumiem pytanie. Cóż, wtedy może po prostu użyć curvelub pointsprzesłonić trzy wykresy;)
nico