Mendapatkan hasil yang berbeda ketika memplot 95% CI elips dengan ggplot atau paket elips

11

Saya ingin memvisualisasikan hasil pengelompokan (diproduksi dengan protoclust{protoclust}) dengan membuat plot scater untuk setiap pasangan variabel yang digunakan untuk mengklasifikasikan data saya, mewarnai dengan kelas dan tumpang tindih elips untuk interval kepercayaan 95% untuk masing-masing kelas (untuk memeriksa yang mana kelas-elip tumpang tindih di bawah setiap pasangan variabel).

Saya telah menerapkan gambar elips dengan dua cara berbeda dan elips yang dihasilkan berbeda! (elips yang lebih besar untuk implementasi pertama!) A priori mereka hanya berbeda dalam ukuran (beberapa skala yang berbeda?), sebagai pusat dan sudut sumbu, tampaknya sama di keduanya. Saya kira saya pasti melakukan sesuatu yang salah dengan menggunakan salah satunya (harap tidak dengan keduanya!), Atau dengan argumen.

Adakah yang bisa memberitahu saya apa yang saya lakukan salah?

Di sini kode untuk dua implementasi; keduanya didasarkan pada jawaban untuk Bagaimana elips data dapat ditumpangkan pada sebar 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)

Berikut adalah dua plot bersama (grafik kiri adalah p1implementasi ( ellipse()):

masukkan deskripsi gambar di sini

Data tersedia di sini: https://www.dropbox.com/sh/xa8xrisa4sfxyj0/l5zaGQmXJt

josetanago
sumber
Ini mungkin tidak masalah, tetapi ketika saya menjalankan kode Anda, saya mendapat peringatan Warning message: In cov.trob(cbind(data$x, data$y)) : Probable convergence failureapakah ini juga terjadi ketika Anda menjalankan kode?
atiretoo - mengembalikan monica
@atiretoo Ya, itu juga terjadi ketika saya menjalankan kode. Saya tidak tahu kenapa. Mungkin orang lain tahu? jose
josetanago

Jawaban:

9

Anda tidak melakukan kesalahan, kedua fungsi tersebut membuat asumsi dasar yang berbeda tentang distribusi data. Implementasi pertama Anda dengan asumsi multivarian normal, dan ke-2 multivariat t-distribusi (lihat? Cov.trob dalam paket MASS). Efeknya lebih mudah dilihat jika Anda menarik satu kelompok:

#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

Jadi meskipun dekat dengan pusat dan orientasi yang sama, mereka tidak sama. Anda dapat mendekati ukuran elips yang sama dengan menggunakan cov.trob()untuk mendapatkan korelasi dan skala untuk lewat ellipse(), dan menggunakan argumen t untuk mengatur skala sama dengan distribusi-f seperti stat_ellipse()halnya.

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

tetapi korespondensi masih belum tepat. Perbedaan harus timbul antara menggunakan dekomposisi cholesky dari matriks kovarians dan menciptakan penskalaan dari korelasi dan standar deviasi. Saya tidak cukup ahli matematika untuk melihat persis di mana perbedaannya.

Yang mana yang benar? Terserah Anda yang memutuskan! The stat_ellipse()implementasi akan kurang sensitif terhadap poin terpencil, sedangkan yang pertama akan lebih konservatif.

atiretoo - mengembalikan monica
sumber
1
Terima kasih banyak telah meluangkan waktu Anda memecahkan pertanyaan ini! Sudah jelas bagiku sekarang. jose
josetanago
1
Saya suka elips di plot, jadi ini menyenangkan untuk melihat fungsi-fungsi ini beraksi.
atiretoo - mengembalikan monica