如何用R/Python手动绘制带VUS值、轴范围0-1的ROC曲面图？

问题背景 #

已知X、Y、Z分别是来自F₁、F₂、F₃累积分布函数（CDF）的随机变量，需手动编写R代码绘制符合要求的ROC曲面图。此前尝试的代码效果不符合预期，原代码如下：

x<-rW(100,7,2)
y<-rW(100,9,4)
z<-rW(100,18,8)
n <- 50
x <- seq(0, 1, length.out = n)
y <- seq(0, 1, length.out = n)
z <- outer(x, y, function(x, y) ifelse(x <= y, x, y))

# Multiplying z values to increase the surface volume
z <- z * 2  # For example, we multiply z values by 2

# Create a 3D surface graph
plot_ly(x = x, y = y, z = z, type = "surface") %>%
  layout(scene = list(
xaxis = list(title = "False Positive Rate (FPR)", range = c(0, 1)),
yaxis = list(title = "True Positive Rate (TPR)", range = c(0, 1)),
zaxis = list(title = "Threshold")
)) 

核心需求：

• ROC曲面图中需显示**VUS（ROC曲面下体积）**值
• 所有坐标轴范围严格限制在0到1之间
• 绘图效果匹配《Statistics in Medicine - 2004 - Nakas - Ordered multiple‐class ROC analysis with continuous measurements》中的图2

正确R实现代码及步骤 #

步骤1：准备环境与生成符合CDF的样本 #

加载必要工具包，生成对应三个类别（F₁、F₂、F₃）的连续型样本，这里用不同参数的正态分布模拟差异CDF：

library(plotly)
library(pracma) # 用于数值积分计算VUS

# 固定随机种子保证结果可复现
set.seed(123)
n_samples <- 200
class1 <- rnorm(n_samples, mean = 1, sd = 0.8)  # 对应F₁
class2 <- rnorm(n_samples, mean = 3, sd = 0.8)  # 对应F₂
class3 <- rnorm(n_samples, mean = 5, sd = 0.8)  # 对应F₃

步骤2：计算ROC曲面的TPR与FPR网格 #

针对有序三分类ROC分析，定义两个阈值t1 ≤ t2，计算对应指标：

• TPR：正确识别类别3的比例（P(X > t2)）
• FPR1：误将类别1判定为2/3的比例（P(X > t1)）
• FPR2：误将类别1/2判定为3的比例（P(X > t2)）

n_grid <- 50
# 生成覆盖所有样本范围的阈值序列
t <- seq(min(c(class1, class2, class3)), max(c(class1, class2, class3)), length.out = n_grid)

# 构建阈值网格，仅保留t1 ≤ t2的有效区域
t1_grid <- rep(t, each = n_grid)
t2_grid <- rep(t, times = n_grid)
mask <- t1_grid <= t2_grid

# 计算各阈值对应的ROC指标
TPR <- sapply(t2_grid[mask], function(t2) mean(class3 > t2))
FPR1 <- sapply(t1_grid[mask], function(t1) mean(class1 > t1))
FPR2 <- sapply(t2_grid[mask], function(t2) mean(class1 > t2 | class2 > t2))

# 将结果整理为矩阵格式，用于曲面绘图
FPR1_mat <- matrix(NA, nrow = n_grid, ncol = n_grid)
FPR2_mat <- matrix(NA, nrow = n_grid, ncol = n_grid)
TPR_mat <- matrix(NA, nrow = n_grid, ncol = n_grid)

FPR1_mat[mask] <- FPR1
FPR2_mat[mask] <- FPR2
TPR_mat[mask] <- TPR

步骤3：计算VUS值 #

通过数值积分计算ROC曲面下的体积：

# 提取有效区域的网格数据
valid_FPR1 <- unique(FPR1)
valid_FPR2 <- unique(FPR2)
valid_TPR <- matrix(TPR, nrow = length(valid_FPR1), ncol = length(valid_FPR2))

# 计算VUS并保留四位小数
VUS <- trapz2d(valid_FPR1, valid_FPR2, valid_TPR)
VUS <- round(VUS, 4)

步骤4：绘制符合要求的ROC曲面图 #

用plotly绘制曲面，锁定坐标轴范围并添加VUS标注：

plot_ly(x = valid_FPR1, y = valid_FPR2, z = valid_TPR, type = "surface",
        colorscale = "Viridis", opacity = 0.8) %>%
  layout(scene = list(
    xaxis = list(title = "FPR (类别1→2/3)", range = c(0, 1)),
    yaxis = list(title = "FPR (类别1/2→3)", range = c(0, 1)),
    zaxis = list(title = "TPR (类别3)", range = c(0, 1)),
    annotations = list(
      x = 0.5, y = 0.5, z = 1,
      text = paste("VUS =", VUS),
      showarrow = FALSE,
      font = list(size = 16, color = "black")
    )
  )) %>%
  add_surface(x = valid_FPR1, y = valid_FPR2, z = matrix(0, nrow = length(valid_FPR1), ncol = length(valid_FPR2)),
              opacity = 0.1, showscale = FALSE) # 添加底面参考，增强立体感

代码说明 #

1. 样本替换：可将正态分布样本替换为真实数据或自定义分布（比如原代码中的rW）
2. ROC逻辑：严格遵循有序三分类ROC的理论定义，仅保留t1 ≤ t2的有效区域，确保曲面形态符合学术规范
3. VUS计算：采用数值积分方法保证精度，结果直观展示在图中
4. 绘图设置：所有坐标轴范围锁定在0-1区间，调整曲面透明度增强可读性，匹配目标文章的图2风格

内容的提问来源于stack exchange，提问作者Ertan Akgenç