Obteniendo resultados diferentes al trazar elipses de IC del 95% con ggplot o el paquete de elipse

11

Quiero visualizar los resultados de una agrupación (producida con protoclust{protoclust}) creando gráficos de dispersión para cada par de variables utilizadas para clasificar mis datos, colorear por clases y superponer las elipses para el intervalo de confianza del 95% para cada una de las clases (para inspeccionar cuáles las clases de elipses se superponen debajo de cada par de variables).

¡He implementado el dibujo de las elipses de dos maneras diferentes y las elipses resultantes son diferentes! (¡elipses más grandes para la primera implementación!) A priori, solo difieren en tamaño (¿alguna escala diferente?), ya que los centros y el ángulo de los ejes parecen ser similares en ambos. Creo que debo estar haciendo algo mal al usar uno de ellos (¡espero que no con ambos!), O con los argumentos.

¿Alguien puede decirme qué estoy haciendo mal?

Aquí el código para las dos implementaciones; ambos se basan en las respuestas a ¿Cómo se puede superponer una elipse de datos en un diagrama de dispersión 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)

Aquí están las dos parcelas juntas (el gráfico izquierdo es p1implementación ( ellipse()):

ingrese la descripción de la imagen aquí

Los datos están disponibles aquí: https://www.dropbox.com/sh/xa8xrisa4sfxyj0/l5zaGQmXJt

josetanago
fuente
Esto podría no importar, pero cuando ejecuto su código recibo una advertencia, Warning message: In cov.trob(cbind(data$x, data$y)) : Probable convergence failure¿esto también sucede cuando ejecuta el código?
atiretoo - reinstalar a monica el
@atiretoo Sí, también sucede cuando ejecuto el código. No se porque. Tal vez alguien más lo sabe? jose
josetanago

Respuestas:

9

No está haciendo nada malo, las dos funciones están haciendo diferentes supuestos subyacentes sobre la distribución de los datos. Su primera implementación supone una multivariada normal y la segunda una distribución t multivariada (consulte? Cov.trob en el paquete MASS). El efecto es más fácil de ver si saca un 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

Entonces, aunque está cerca del mismo centro y orientación, no son lo mismo. Puede acercarse a la elipse del mismo tamaño utilizando cov.trob()para obtener la correlación y la escala para pasar ellipse(), y utilizando el argumento t para establecer la escala igual a una distribución f como lo stat_ellipse()hace.

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

pero la correspondencia aún no es exacta. La diferencia debe surgir entre el uso de la descomposición cholesky de la matriz de covarianza y la creación de la escala de la correlación y las desviaciones estándar. No soy lo suficientemente matemático para ver exactamente dónde está la diferencia.

¿Cuál es el correcto? ¡Eso depende de ti decidir! La stat_ellipse()implementación será menos sensible a los puntos periféricos, mientras que la primera será más conservadora.

atiretoo - restablecer monica
fuente
1
¡Muchas gracias por tomarse su tiempo para resolver esta pregunta! Está claro para mí ahora. jose
josetanago
1
Me gustan las elipses en las tramas, así que fue divertido ver estas funciones en acción.
atiretoo - reinstalar a monica el