Uzyskanie różnych wyników podczas kreślenia elips 95% CI za pomocą ggplot lub pakietu elipsy

11

Chcę wizualizować wyniki grupowania (wytworzone za pomocą protoclust{protoclust}) poprzez tworzenie wykresów scater dla każdej pary zmiennych używanych do klasyfikowania moich danych, kolorowania według klas i nakładania się elips na 95% przedział ufności dla każdej z klas (aby sprawdzić, które klasy elipses pokrywają się pod każdą parą zmiennych).

Zaimplementowałem rysunek elips na dwa różne sposoby, a powstałe elipsy są różne! (większe elipsy dla pierwszej implementacji!) A priori różnią się jedynie rozmiarem (niektóre różne skalowanie?), ponieważ środki i kąt osi wydają się być podobne w obu. Chyba muszę coś zrobić źle, używając jednego z nich (mam nadzieję, że nie w obu przypadkach!) Lub w argumentach.

Czy ktoś może mi powiedzieć, co robię źle?

Tutaj kod dla dwóch implementacji; oba są oparte na odpowiedziach na pytanie, w jaki sposób elipsę danych można nałożyć na wykres punktowy ggplot2?

### 1st implementation 
### using ellipse{ellipse}
library(ellipse)
library(ggplot2) 
library(RColorBrewer)
colorpal <- brewer.pal(10, "Paired")

x <- data$x
y <- data$y
group <- data$group
df <- data.frame(x=x, y=y, group=factor(group))

df_ell <- data.frame() 
for(g in levels(df$group)){df_ell <- rbind(df_ell, cbind(as.data.frame(with(df[df$group==g,], ellipse(cor(x, y),scale=c(sd(x),sd(y)),centre=c(mean(x),mean(y))))),group=g))} 

p1 <- ggplot(data=df, aes(x=x, y=y,colour=group)) + geom_point() + 
  geom_path(data=df_ell, aes(x=x, y=y,colour=group))+scale_colour_manual(values=colorpal)

### 2nd implementation 
###using function ellipse_stat() 
###code by Josef Fruehwald available in: https://github.com/JoFrhwld/FAAV/blob/master/r/stat-ellipse.R

p2 <-qplot(data=df, x=x,y=y,colour=group)+stat_ellipse(level=0.95)+scale_colour_manual(values=colorpal)

Oto dwa wykresy razem (lewy wykres to p1implementacja ( ellipse()):

wprowadź opis zdjęcia tutaj

Dane są dostępne tutaj: https://www.dropbox.com/sh/xa8xrisa4sfxyj0/l5zaGQmXJt

josetanago
źródło
To może nie mieć znaczenia, ale kiedy uruchamiam kod, pojawia się ostrzeżenie, Warning message: In cov.trob(cbind(data$x, data$y)) : Probable convergence failureczy dzieje się tak również wtedy, gdy kod jest uruchamiany?
atiretoo
@atiretoo Tak, dzieje się tak również po uruchomieniu kodu. Nie wiem dlaczego. Może ktoś jeszcze wie? jose
josetanago

Odpowiedzi:

9

Nie robisz nic złego, obie funkcje przyjmują różne podstawowe założenia dotyczące dystrybucji danych. Twoja pierwsza implementacja zakłada wielowymiarową normalność, a druga wielowymiarową dystrybucję t (patrz? Cov.trob w pakiecie MASS). Efekt łatwiej zobaczyć, jeśli wyciągniesz jedną grupę:

#pull out group 1
pick = group ==1
p3 <- qplot(data=df[pick,], x=x, y=y)
tl = with(df[pick,], 
     ellipse(cor(x, y),scale=c(sd(x),sd(y)),
             centre=c(mean(x),mean(y))))
p3 <- p3 + geom_path(data=as.data.frame(tl), aes(x=x, y=y))
p3 <- p3 + stat_ellipse(level=0.95)
p3 # looks off center
p3 <- p3 + geom_point(aes(x=mean(x),y=mean(y),size=2,color="red"))
p3

Więc chociaż jest blisko tego samego centrum i orientacji, nie są one takie same. Można zbliżyć się do tego samego rozmiaru za pomocą elipsy cov.trob(), aby uzyskać korelację i skalę dla przechodząc do ellipse(), i używając argumentu t ustawić skalowanie równa f-dystrybucji, jak stat_ellipse()robi.

tcv = cov.trob(data[pick,2:3],cor=TRUE)
tl = with(df[pick,], 
          ellipse(tcv$cor[2,1],scale=sqrt(diag(tcv$cov)),
                  t=qf(0.95,2,length(x)-1),
                  centre=tcv$center))
p3 <- p3 + geom_path(data=as.data.frame(tl), aes(x=x, y=y,color="red"))
p3

ale korespondencja wciąż nie jest dokładna. Różnica musi wynikać z zastosowania chłodnego rozkładu macierzy kowariancji i utworzenia skalowania na podstawie korelacji i odchyleń standardowych. Nie jestem wystarczająco matematykiem, żeby zobaczyć, gdzie dokładnie jest różnica.

Który jest prawidłowy? Ty decydujesz! stat_ellipse()Realizacja będzie mniej wrażliwy na odległych punktów, podczas gdy pierwszy będzie bardziej konserwatywna.

atiretoo - przywróć monikę
źródło
1
Wielkie dzięki za poświęcenie czasu na rozwiązanie tego pytania! Teraz jest dla mnie jasne. jose
josetanago
1
Lubię elipsy na wykresach, więc fajnie było zobaczyć te funkcje w akcji.
atiretoo