paso - R: trazado de probabilidades de clasificación posterior de un análisis discriminante lineal en ggplot2

analisis discriminante paso a paso

Question

Usando ggord se puede hacer un buen análisis discriminante lineal ggplot2 biplots (ver capítulo 11, Fig. 11.5 en "Biplots en la práctica" por M. Greenacre), como en

library(MASS)
install.packages("devtools")
library(devtools)
install_github("fawda123/ggord")
library(ggord)
data(iris)
ord <- lda(Species ~ ., iris, prior = rep(1, 3)/3)
ggord(ord, iris$Species)

También me gustaría agregar las regiones de clasificación (mostradas como regiones sólidas del mismo color que su grupo respectivo con, por ejemplo, alfa = 0.5) o las probabilidades posteriores de pertenencia a una clase (con alfa que varía según la probabilidad posterior y el mismo color que utilizado para cada grupo) (como se puede hacer en BiplotGUI , pero estoy buscando una solución ggplot2 ). ¿Alguien sabría cómo hacer esto con ggplot2 , quizás usando geom_tile ?

EDITAR: debajo, alguien pregunta cómo calcular las probabilidades de clasificación posterior y las clases predichas. Esto es así:

library(MASS)
library(ggplot2)
library(scales)
fit <- lda(Species ~ ., data = iris, prior = rep(1, 3)/3)
datPred <- data.frame(Species=predict(fit)$class,predict(fit)$x)
#Create decision boundaries
fit2 <- lda(Species ~ LD1 + LD2, data=datPred, prior = rep(1, 3)/3)
ld1lim <- expand_range(c(min(datPred$LD1),max(datPred$LD1)),mul=0.05)
ld2lim <- expand_range(c(min(datPred$LD2),max(datPred$LD2)),mul=0.05)
ld1 <- seq(ld1lim[[1]], ld1lim[[2]], length.out=300)
ld2 <- seq(ld2lim[[1]], ld1lim[[2]], length.out=300)
newdat <- expand.grid(list(LD1=ld1,LD2=ld2))
preds <-predict(fit2,newdata=newdat)
predclass <- preds$class
postprob <- preds$posterior
df <- data.frame(x=newdat$LD1, y=newdat$LD2, class=predclass)
df$classnum <- as.numeric(df$class)
df <- cbind(df,postprob)
head(df)

           x        y     class classnum       setosa   versicolor virginica
1 -10.122541 -2.91246 virginica        3 5.417906e-66 1.805470e-10         1
2 -10.052563 -2.91246 virginica        3 1.428691e-65 2.418658e-10         1
3  -9.982585 -2.91246 virginica        3 3.767428e-65 3.240102e-10         1
4  -9.912606 -2.91246 virginica        3 9.934630e-65 4.340531e-10         1
5  -9.842628 -2.91246 virginica        3 2.619741e-64 5.814697e-10         1
6  -9.772650 -2.91246 virginica        3 6.908204e-64 7.789531e-10         1

colorfun <- function(n,l=65,c=100) { hues = seq(15, 375, length=n+1); hcl(h=hues, l=l, c=c)[1:n] } # default ggplot2 colours
colors <- colorfun(3)
colorslight <- colorfun(3,l=90,c=50)
ggplot(datPred, aes(x=LD1, y=LD2) ) +
    geom_raster(data=df, aes(x=x, y=y, fill = factor(class)),alpha=0.7,show_guide=FALSE) +
    geom_contour(data=df, aes(x=x, y=y, z=classnum), colour="red2", alpha=0.5, breaks=c(1.5,2.5)) +
    geom_point(data = datPred, size = 3, aes(pch = Species,  colour=Species)) +
    scale_x_continuous(limits = ld1lim, expand=c(0,0)) +
    scale_y_continuous(limits = ld2lim, expand=c(0,0)) +
    scale_fill_manual(values=colorslight,guide=F)

(no estoy totalmente seguro de que este enfoque para mostrar los límites de clasificación utilizando contornos / roturas en 1.5 y 2.5 sea siempre correcto; es correcto para el límite entre las especies 1 y 2 y las especies 2 y 3, pero no si la región de la especie 1 sería al lado de la especie 3, ya que obtendría dos límites allí entonces, tal vez tendría que usar el enfoque utilizado here donde cada límite entre cada par de especies se considera por separado)

Esto me lleva a trazar las regiones de clasificación. Sin embargo, estoy buscando una solución para trazar las probabilidades reales de clasificación posterior para cada especie en cada coordenada, utilizando alfa (opacidad) proporcional a la probabilidad de clasificación posterior para cada especie y un color específico de la especie. En otras palabras, con una pila de tres imágenes superpuestas. Como se sabe que la mezcla alfa en ggplot2 order-dependent , creo que los colores de esta pila tendrían que calcularse de antemano, sin embargo, y graficar usando algo como

qplot(x, y, data=mydata, fill=rgb, geom="raster") + scale_fill_identity() 

Aquí hay un ejemplo de SAS de lo que busco :

¿Alguien sabría cómo hacer esto tal vez? ¿O alguien tiene alguna idea sobre cómo representar mejor estas probabilidades de clasificación posterior?

Tenga en cuenta que el método debería funcionar para cualquier número de grupos, no solo para este ejemplo específico.

Answer 1

Supongo que la forma más fácil será mostrar las probabilidades posteriores. Es bastante sencillo para su caso:

datPred$maxProb <- apply(predict(fit)$posterior, 1, max)
ggplot(datPred, aes(x=LD1, y=LD2) ) +
  geom_raster(data=df, aes(x=x, y=y, fill = factor(class)),alpha=0.7,show_guide=FALSE) +
  geom_contour(data=df, aes(x=x, y=y, z=classnum), colour="red2", alpha=0.5, breaks=c(1.5,2.5)) +
  geom_point(data = datPred, size = 3, aes(pch = Species,  colour=Species, alpha = maxProb)) +
  scale_x_continuous(limits = ld1lim, expand=c(0,0)) +
  scale_y_continuous(limits = ld2lim, expand=c(0,0)) +
  scale_fill_manual(values=colorslight, guide=F)

Puedes ver los puntos mezclados en el borde azul-verde.

Answer 2

También se le ocurrió la siguiente solución fácil: simplemente haga una columna en df donde las predicciones de clase se realicen de forma estocástica, de acuerdo con las probabilidades posteriores, lo que resultará en interpolación en regiones inciertas, por ejemplo, como en

fit = lda(Species ~ Sepal.Length + Sepal.Width, data = iris, prior = rep(1, 3)/3)
ld1lim <- expand_range(c(min(datPred$LD1),max(datPred$LD1)),mul=0.5)
ld2lim <- expand_range(c(min(datPred$LD2),max(datPred$LD2)),mul=0.5)

descansa como arriba, e insertando

lvls=unique(df$class)
df$classpprob=apply(df[,as.character(lvls)],1,function(row) sample(lvls,1,prob=row))

p=ggplot(datPred, aes(x=LD1, y=LD2) ) +
  geom_raster(data=df, aes(x=x, y=y, fill = factor(classpprob)),hpad=0, vpad=0, alpha=0.7,show_guide=FALSE) +
  geom_point(data = datPred, size = 3, aes(pch = Group,  colour=Group)) +
  scale_fill_manual(values=colorslight,guide=F) +
  scale_x_continuous(limits=rngs[[1]], expand=c(0,0)) +
  scale_y_continuous(limits=rngs[[2]], expand=c(0,0))

me da

Mucho más fácil y claro que comenzar a mezclar colores de alguna manera aditiva o sustractiva de todos modos (que es la parte en la que todavía tenía problemas, y que aparentemente no es tan trivial hacerlo bien).