Obtendo resultados diferentes ao plotar elipses de IC de 95% com o ggplot ou o pacote ellipse

11

Quero visualizar os resultados de um agrupamento (produzido com protoclust{protoclust}) criando gráficos de dispersão para cada par de variáveis ​​usadas para classificar meus dados, colorindo por classes e sobrepondo as elipses para o intervalo de confiança de 95% para cada uma das classes (para inspecionar quais classes de elipses se sobrepõem sob cada par de variáveis).

Eu implementei o desenho das elipses de duas maneiras diferentes e as elipses resultantes são diferentes! (elipses maiores para a primeira implementação!) A priori, diferem apenas em tamanho (algumas escalas diferentes?), pois o centro e o ângulo dos eixos parecem os mesmos em ambos. Acho que devo estar fazendo algo errado usando um deles (espero que não com os dois!), Ou com os argumentos.

Alguém pode me dizer o que estou fazendo de errado?

Aqui o código para as duas implementações; ambos são baseados nas respostas de Como uma elipse de dados pode ser sobreposta em um gráfico de dispersão 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)

Aqui estão os dois gráficos juntos (o gráfico à esquerda é p1implementação ( ellipse()):

insira a descrição da imagem aqui

Os dados estão disponíveis aqui: https://www.dropbox.com/sh/xa8xrisa4sfxyj0/l5zaGQmXJt

josetanago
fonte
Isso pode não importar, mas quando executo seu código, recebo um aviso. Warning message: In cov.trob(cbind(data$x, data$y)) : Probable convergence failureIsso também acontece quando você executa o código?
atiretoo - reinstate monica
@atiretoo Sim, isso também acontece quando executo o código. Eu não sei porque. Talvez alguém mais saiba? jose
josetanago

Respostas:

9

Você não está fazendo nada errado, as duas funções estão fazendo suposições subjacentes diferentes sobre a distribuição dos dados. Sua primeira implementação está assumindo multivariada normal e a segunda uma distribuição t multivariada (consulte? Cov.trob no pacote MASS). O efeito é mais fácil de ver se você seleciona um grupo:

#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

Portanto, embora esteja próximo do mesmo centro e orientação, eles não são os mesmos. Você pode se aproximar da elipse do mesmo tamanho usando cov.trob()para obter a correlação e a escala para passar para ellipse()e usando o argumento t para definir o dimensionamento igual a uma distribuição f como stat_ellipse()faz.

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

mas a correspondência ainda não é exata. A diferença deve estar surgindo entre o uso da decomposição de Cholesky da matriz de covariância e a criação do dimensionamento a partir da correlação e dos desvios padrão. Não sou matemático o suficiente para ver exatamente onde está a diferença.

Qual deles está correto? Cabe a você decidir! A stat_ellipse()implementação será menos sensível aos pontos periféricos, enquanto a primeira será mais conservadora.

atiretoo - restabelecer monica
fonte
1
Muito obrigado por dedicar seu tempo resolvendo esta pergunta! Está claro para mim agora. jose
josetanago
1
Eu gosto de elipses em gráficos, então foi divertido ver essas funções em ação.
atiretoo - restabelece monica