如何在R的3D图中绘制分类算法生成的分界面

我尝试在R中使用plot3D包绘制分类算法得到的分区边界。如果仅有两个预测变量，仅需两个坐标轴即可完成绘制（例如使用partimat函数），任务相对简单，但我至今未找到基于三个预测变量的3D空间分类分界面的理想绘制方案。

2D场景参考实现 #

为直观说明问题，我们先在鸢尾花数据集上使用线性判别分析（LDA）算法构建仅包含两个轴的分区：

# 加载包并裁剪鸢尾花数据集
library(klaR)

data = droplevels(iris[iris$Species != 'virginica', ])

partimat(Species ~ Sepal.Length + Sepal.Width, data, 
         method = 'lda')

运行后会得到一张2D图，两个鸢尾花品种之间有清晰的分区边界：

3D场景现有实现及问题 #

但partimat一次最多只能处理两个预测变量（可参考?partimat文档）。接下来我们看3D场景的问题：

library(plot3D)
    
# 绘制原始数据
points3D(data$Sepal.Length, data$Sepal.Width, data$Petal.Length,
             colkey = F,
             pch = 16, cex = 2,
             theta = 30, phi = 30, 
             ticktype = 'detailed',
             col = data$Species)

我希望基于LDA这类分类算法绘制分隔两类数据的平面。参考现有示例，我尝试了一种实现方案：先构建包含三个预测变量的LDA模型，再在三个预测变量的最大、最小值范围内生成一系列网格点，预测这些点对应的品种（山鸢尾或变色鸢尾），将预测结果以不同颜色绘制在3D图中，就能得到表示不同品种对应3D空间的点云，分区可以通过点云中红蓝点的模糊交界看出来：

# 构建包含三个预测变量的分类模型
m = lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length, data)

# 对三个特征的全量取值范围预测分类标签
np = 50

nx = seq(from = min(data[, 1]), to = max(data[, 1]), length.out = np)
ny = seq(from = min(data[, 2]), to = max(data[, 2]), length.out = np)
nz = seq(from = min(data[, 3]), to = max(data[, 3]), length.out = np)
nd = expand.grid(Sepal.Length = nx, Sepal.Width = ny, Petal.Length = nz)

p    = as.numeric(predict(m, newdata = nd)$class)
part = cbind(nd, Partition = p)

# 绘制分区点云并叠加原始数据点
scatter3D(part$Sepal.Length, part$Sepal.Width, part$Petal.Length, 
          colvar = part$Partition, 
          colkey = F,
          alpha = 0.5,
          pch = 16, cex = 0.3, 
          theta = 30, phi = 30, 
          ticktype = 'detailed',
          plot = F)
points3D(data$Sepal.Length, data$Sepal.Width, data$Petal.Length,
         colkey = F,
         pch = 16, cex = 2,
         theta = 30, phi = 30, 
         ticktype = 'detailed',
         col = data$Species,
         add = T)

但这个方案并不理想：点云会遮挡数据点，且本身也比较干扰视觉。即使调整点的透明度也只能略微改善效果，我认为更好的方案是在两类的交界位置绘制一个平面（类似回归平面）。另外我最终还需要适配随机森林等不同分类器，不希望方案仅支持LDA这类线性模型。

通用方案：适配所有分类器的等值面绘制 #

使用misc3d包的contour3d函数提取分类分界的等值面，直接绘制半透明分界面，完全避免点云的遮挡问题。该方案不限制分类器类型，LDA、随机森林、SVM等所有能输出分类概率的模型都可以直接适配。

实现代码 #

# 加载依赖包
library(klaR)
library(MASS)
library(plot3D)
library(misc3d)

# 数据准备与模型训练（可替换为任意分类模型）
data = droplevels(iris[iris$Species != 'virginica', ])
# 此处可替换为randomForest、svm等模型训练代码
m = lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length, data)

# 生成三维网格
np = 50
nx = seq(from = min(data[, 1]), to = max(data[, 1]), length.out = np)
ny = seq(from = min(data[, 2]), to = max(data[, 2]), length.out = np)
nz = seq(from = min(data[, 3]), to = max(data[, 3]), length.out = np)
nd = expand.grid(Sepal.Length = nx, Sepal.Width = ny, Petal.Length = nz)

# 计算网格点属于第二类的概率，生成分界阈值三维数组
nd$prob = predict(m, newdata = nd)$posterior[,2]
prob_arr = array(nd$prob, dim = c(np, np, np))

# 绘制原始数据点
points3D(data$Sepal.Length, data$Sepal.Width, data$Petal.Length,
             colkey = F,
             pch = 16, cex = 2,
             theta = 30, phi = 30, 
             ticktype = 'detailed',
             col = data$Species)

# 叠加半透明分类分界面（阈值设为0.5，即两类分界）
contour3d(x = nx, y = ny, z = nz, f = prob_arr, level = 0.5, 
          alpha = 0.3, col = "gray", add = TRUE)

方案优势 #

• 半透明分界面不会遮挡原始数据点，视觉效果清爽
• 更换分类器仅需修改模型训练和概率预测的代码，其余逻辑完全复用
• 支持多分类场景，设置多个level参数即可绘制多个分界面

线性分类器专属简化方案 #

如果仅使用LDA、逻辑回归这类线性分类器，可以直接推导决策面的线性方程，用planes3D函数直接绘制平面，运算效率更高。

# 提取LDA模型的决策面系数
coef = m$scaling
intercept = -mean(m$means %*% coef)

# 绘制原始数据点
points3D(data$Sepal.Length, data$Sepal.Width, data$Petal.Length,
             colkey = F,
             pch = 16, cex = 2,
             theta = 30, phi = 30, 
             ticktype = 'detailed',
             col = data$Species)

# 直接绘制线性分界面
planes3D(a = coef[1], b = coef[2], c = coef[3], d = intercept,
         alpha = 0.3, col = "gray", add = TRUE)

内容的提问来源于stack exchange，提问作者Cai Ladd