R：在 ggplot2 中绘制线性判别分析的后验分类概率
R: plotting posterior classification probabilities of a linear discriminant analysis in ggplot2
使用 ggord 可以进行很好的线性判别分析 ggplot2 双标图（参见 M. Greenacre "Biplots in practice" 中的第 11 章图 11.5），如
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)
我还想添加 class 化区域（显示为与其各自组颜色相同的实心区域，例如 alpha = 0.5）或 class 成员资格的后验概率（ alpha 然后根据这个后验概率和每组使用的相同颜色而变化）（可以在 BiplotGUI 中完成，但我正在寻找 ggplot2 解决方案）。有谁知道如何使用 ggplot2，也许使用 geom_tile？
编辑：下面有人问如何计算后验class化概率和预测classes。事情是这样的：
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)
（不太确定这种使用 contours/breaks 在 1.5 和 2.5 处显示 class 化边界的方法总是正确的 - 它对于物种 1 和 2 以及物种 2 和3，但如果物种 1 的区域紧挨着物种 3，则不会，因为那时我会在那里得到两个边界 - 也许我将不得不使用使用的方法 here，其中每个物种对之间的每个边界都被单独考虑)
这让我可以绘制 class化区域。我正在寻找一种解决方案，尽管它还可以绘制每个物种在每个坐标处的实际后验 class 化概率，使用与每个物种的后验 class 化概率成比例的 alpha（不透明性），以及一个物种-特定颜色。换句话说，叠加了三个图像。由于 ggplot2 中的 alpha 混合已知为 order-dependent，我认为此堆栈的颜色必须事先计算，并使用类似
的方法绘制
qplot(x, y, data=mydata, fill=rgb, geom="raster") + scale_fill_identity()
Here is a SAS example of what I am after:
有人知道怎么做吗？或者有人对如何最好地表示这些后验 class 化概率有任何想法吗？
请注意，该方法适用于任意数量的组，而不仅仅是这个特定示例。
我想最简单的方法是显示后验概率。这对你的情况来说非常简单：
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)
您可以看到点在蓝绿色边框处混合。
还提出了以下简单的解决方案：只需在 df 中创建一个列，其中 class 预测是根据后验概率随机进行的，然后会导致不确定区域的抖动，例如如
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)
如上休息，插入
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))
给我
无论如何，比开始以某种加色或减色方式混合颜色要容易和清晰得多（这是我仍然遇到问题的部分，而且显然要做好也不是那么微不足道）。
使用 ggord 可以进行很好的线性判别分析 ggplot2 双标图（参见 M. Greenacre "Biplots in practice" 中的第 11 章图 11.5），如
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)
我还想添加 class 化区域（显示为与其各自组颜色相同的实心区域，例如 alpha = 0.5）或 class 成员资格的后验概率（ alpha 然后根据这个后验概率和每组使用的相同颜色而变化）（可以在 BiplotGUI 中完成，但我正在寻找 ggplot2 解决方案）。有谁知道如何使用 ggplot2，也许使用 geom_tile？
编辑：下面有人问如何计算后验class化概率和预测classes。事情是这样的：
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)
（不太确定这种使用 contours/breaks 在 1.5 和 2.5 处显示 class 化边界的方法总是正确的 - 它对于物种 1 和 2 以及物种 2 和3，但如果物种 1 的区域紧挨着物种 3，则不会，因为那时我会在那里得到两个边界 - 也许我将不得不使用使用的方法 here，其中每个物种对之间的每个边界都被单独考虑)
这让我可以绘制 class化区域。我正在寻找一种解决方案，尽管它还可以绘制每个物种在每个坐标处的实际后验 class 化概率，使用与每个物种的后验 class 化概率成比例的 alpha（不透明性），以及一个物种-特定颜色。换句话说，叠加了三个图像。由于 ggplot2 中的 alpha 混合已知为 order-dependent，我认为此堆栈的颜色必须事先计算，并使用类似
的方法绘制qplot(x, y, data=mydata, fill=rgb, geom="raster") + scale_fill_identity()
Here is a SAS example of what I am after:
有人知道怎么做吗？或者有人对如何最好地表示这些后验 class 化概率有任何想法吗？
请注意，该方法适用于任意数量的组，而不仅仅是这个特定示例。
我想最简单的方法是显示后验概率。这对你的情况来说非常简单：
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)
您可以看到点在蓝绿色边框处混合。
还提出了以下简单的解决方案：只需在 df 中创建一个列，其中 class 预测是根据后验概率随机进行的，然后会导致不确定区域的抖动，例如如
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)
如上休息，插入
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))
给我
无论如何，比开始以某种加色或减色方式混合颜色要容易和清晰得多（这是我仍然遇到问题的部分，而且显然要做好也不是那么微不足道）。