我想在我的三维依赖图中，根据拟合值添加一系列渐变颜色（例如，较高的拟合值使用较深的颜色，较低的拟合值使用较浅的颜色）。
我使用了dismo包中提供的示例：
library(dismo) data(Anguilla_train)angaus.tc5.lr01 <- gbm.step(data=Anguilla_train, gbm.x = 3:13, gbm.y = 2,family = "bernoulli", tree.complexity = 5, learning.rate = 0.01, bag.fraction = 0.5)# 在gbm模型中查找交互作用：find.int <- gbm.interactions( angaus.tc5.lr01)find.int$interactionsfind.int$rank.list我只能将整个图设置为同一种颜色：
gbm.perspec( angaus.tc5.lr01, 7, 1, x.label = "USRainDays", y.label = "SegSumT", z.label = "Fitted values", z.range=c(0,0.435), col="blue")
或者添加渐变色，但不依据拟合值：
gbm.perspec( angaus.tc5.lr01, 7, 1, x.label = "USRainDays", y.label = "SegSumT", z.label = "Fitted values", col=heat.colors(50), z.range=c(0,0.435))
我还检查了gbm.perspec函数的代码，如果我理解正确的话，拟合值在公式中被称为“prediction”，后来成为了传递给最终绘图的“pred.matrix”的一部分：persp(x = x.var, y = y.var, z = pred.matrix…)，但我无法从gbm.perspec公式中访问它们。我尝试通过在函数内的persp()中添加“col=heat.colors(100)[round(pred.matrix*100, 0)]”来修改gbm.perpec函数，但这并未达到我想要的效果：
persp(x = x.var, y = y.var, z = pred.matrix, zlim = z.range, xlab = x.label, ylab = y.label, zlab = z.label, theta = theta, phi = phi, r = sqrt(10), d = 3, ticktype = ticktype, col=heat.colors(100)[round(pred.matrix*100, 0)], mgp = c(4, 1, 0), ...)
我认为解决方案可能来自于修改gbm.perpec函数，您知道如何做吗？
感谢您的宝贵时间！
回答：
通过将这里找到的一些代码整合到gbm.perspec()的源代码中，您可以创建所需的效果。
首先运行
# 颜色调色板（100种颜色）col.pal<-colorRampPalette(c("blue", "red"))colors<-col.pal(100)然后，在gbm.perspec()源代码中的else之后添加z.facet.center，并将代码中的z更改为pred.matrix，如下所示，
# 最后绘制结果#if (!perspective) { image(x = x.var, y = y.var, z = pred.matrix, zlim = z.range)} else { z.facet.center <- (pred.matrix[-1, -1] + pred.matrix[-1, -ncol(pred.matrix)] + pred.matrix[-nrow(pred.matrix), -1] + pred.matrix[-nrow(pred.matrix), -ncol(pred.matrix)])/4 # 面中心的范围按100比例（颜色数量） z.facet.range<-cut(z.facet.center, 100) persp(x=x.var, y=y.var, z=pred.matrix, zlim= z.range, # 输入变量 xlab = x.label, ylab = y.label, zlab = z.label, # 标签 theta=theta, phi=phi, r = sqrt(10), d = 3, col=colors[z.facet.range],# 视图参数 ticktype = ticktype, mgp = c(4,1,0), ...) #这将为您生成如下图表（请注意，这不是使用示例数据集绘制的，因此交互效果与问题中的图表不同）。 
或者，您可以创建一个新函数。以下示例修改了gbm.perspec()以提供从白色到红色的渐变。只需在R中运行代码，然后将gbm.perspec()更改为gbm.perspec2()
# 交互函数 # 颜色调色板（100种颜色）col.pal<-colorRampPalette(c("white", "pink", "red"))colors<-col.pal(100)gbm.perspec2 <- function(gbm.object, x = 1, # 要绘制的第一个变量 y = 2, # 要绘制的第二个变量 pred.means = NULL, # 允许为其他变量指定值 x.label = NULL, # 允许手动指定x轴标签 x.range = NULL, # 手动指定x变量的范围 y.label = NULL, # 以及y轴标签 z.label = "fitted value", #默认z轴标签 y.range = NULL, # 以及y z.range = NULL, # 允许控制垂直轴 leg.coords = NULL, #可以指定图例的坐标（x, y） ticktype = "detailed",# 指定详细类型 - 否则为"simple" theta = 55, # 旋转角度 phi=40, # 和仰角 smooth = "none", # 控制预测表面的平滑 mask = FALSE, # 控制是否使用样本强度模型进行掩蔽 perspective = TRUE, # 控制是绘制轮廓图还是透视图 ...) # 允许传递额外的参数给绘图例程# 有用的选项包括shade, ltheta, lphi来控制照明# 和cex来控制文本大小 - cex.axis和cex.lab没有效果{ if (! requireNamespace('gbm') ) { stop('您需要安装gbm包来使用此函数') } requireNamespace('splines') #获取提升模型的详细信息 gbm.call <- gbm.object$gbm.call gbm.x <- gbm.call$gbm.x n.preds <- length(gbm.x) gbm.y <- gbm.call$gbm.y pred.names <- gbm.call$predictor.names family = gbm.call$family # 现在为x和y预测变量设置范围变量 have.factor <- FALSE x.name <- gbm.call$predictor.names[x] if (is.null(x.label)) { x.label <- gbm.call$predictor.names[x] } y.name <- gbm.call$predictor.names[y] if (is.null(y.label)) { y.label <- gbm.call$predictor.names[y] } data <- gbm.call$dataframe[ , gbm.x, drop=FALSE] n.trees <- gbm.call$best.trees # 如果边缘变量是向量，则在范围内创建间隔 if (is.vector(data[,x])) { if (is.null(x.range)) { x.var <- seq(min(data[,x],na.rm=T),max(data[,x],na.rm=T),length = 50) } else { x.var <- seq(x.range[1],x.range[2],length = 50) } } else { x.var <- names(table(data[,x])) have.factor <- TRUE } if (is.vector(data[,y])) { if (is.null(y.range)) { y.var <- seq(min(data[,y],na.rm=T),max(data[,y],na.rm=T),length = 50) } else {y.var <- seq(y.range[1],y.range[2],length = 50)} } else { y.var <- names(table(data[,y])) if (have.factor) { #检查我们是否已经有一个因子 stop("至少有一个边缘预测变量必须是向量！") } else {have.factor <- TRUE} } pred.frame <- expand.grid(list(x.var,y.var)) names(pred.frame) <- c(x.name,y.name) pred.rows <- nrow(pred.frame) #确保因子变量首先出现 if (have.factor) { if (is.factor(pred.frame[,2])) { # 交换它们 pred.frame <- pred.frame[,c(2,1)] x.var <- y.var } } j <- 3 # 循环遍历预测变量 # 如果是非目标变量，找到均值 for (i in 1:n.preds) { if (i != x & i != y) { if (is.vector(data[,i])) { m <- match(pred.names[i],names(pred.means)) if (is.na(m)) { pred.frame[,j] <- mean(data[,i],na.rm=T) } else pred.frame[,j] <- pred.means[m] } if (is.factor(data[,i])) { m <- match(pred.names[i],names(pred.means)) temp.table <- table(data[,i]) if (is.na(m)) { pred.frame[,j] <- rep(names(temp.table)[2],pred.rows) } else { pred.frame[,j] <- pred.means[m] } pred.frame[,j] <- factor(pred.frame[,j],levels=names(temp.table)) } names(pred.frame)[j] <- pred.names[i] j <- j + 1 } } # # 形成预测 # #assign("pred.frame", pred.frame, pos=1) prediction <- gbm::predict.gbm(gbm.object,pred.frame,n.trees = n.trees, type="response") #assign("prediction", prediction, pos=1, immediate =T) # 如果需要，模型平滑 if (smooth == "model") { pred.glm <- glm(prediction ~ ns(pred.frame[,1], df = 8) * ns(pred.frame[,2], df = 8), data=pred.frame,family=poisson) prediction <- fitted(pred.glm) } # 报告最大值并为z设置现实范围 max.pred <- max(prediction) message("最大值 = ",round(max.pred,2),"\n") if (is.null(z.range)) { if (family == "bernoulli") { z.range <- c(0,1) } else if (family == "poisson") { z.range <- c(0,max.pred * 1.1) } else { z.min <- min(data[,y],na.rm=T) z.max <- max(data[,y],na.rm=T) z.delta <- z.max - z.min z.range <- c(z.min - (1.1 * z.delta), z.max + (1.1 * z.delta)) } } # 现在处理假设x和y都是向量 if (have.factor == FALSE) { # 形成矩阵 pred.matrix <- matrix(prediction,ncol=50,nrow=50) # 如果需要，进行核平滑 if (smooth == "average") { #应用3x3平滑平均 pred.matrix.smooth <- pred.matrix for (i in 2:49) { for (j in 2:49) { pred.matrix.smooth[i,j] <- mean(pred.matrix[c((i-1):(i+1)),c((j-1):(j+1))]) } } pred.matrix <- pred.matrix.smooth } # 屏蔽超矩形内但样本空间外的值 if (mask) { mask.trees <- gbm.object$gbm.call$best.trees point.prob <- gbm::predict.gbm(gbm.object[[1]],pred.frame, n.trees = mask.trees, type="response") point.prob <- matrix(point.prob,ncol=50,nrow=50) pred.matrix[point.prob < 0.5] <- 0.0 } # # 最后绘制结果 # if (!perspective) { image(x = x.var, y = y.var, z = pred.matrix, zlim = z.range) } else { z.facet.center <- (pred.matrix[-1, -1] + pred.matrix[-1, -ncol(pred.matrix)] + pred.matrix[-nrow(pred.matrix), -1] + pred.matrix[-nrow(pred.matrix), -ncol(pred.matrix)])/4 # 面中心的范围按100比例（颜色数量） z.facet.range<-cut(z.facet.center, 100) persp(x=x.var, y=y.var, z=pred.matrix, zlim= z.range, # 输入变量 xlab = x.label, ylab = y.label, zlab = z.label, # 标签 theta=theta, phi=phi, r = sqrt(10), d = 3, col=colors[z.facet.range],# 视图参数 ticktype = ticktype, mgp = c(4,1,0), ...) # } } if (have.factor) { # 我们需要为每个x绘制y的值 factor.list <- names(table(pred.frame[,1])) n <- 1 #添加这一部分，以便z.range仍然按预期工作： if (is.null(z.range)) { vert.limits <- c(0, max.pred * 1.1) } else { vert.limits <- z.range } plot(pred.frame[pred.frame[,1]==factor.list[1],2], prediction[pred.frame[,1]==factor.list[1]], type = 'l', #ylim = c(0, max.pred * 1.1), ylim = vert.limits, xlab = y.label, ylab = z.label, ...) for (i in 2:length(factor.list)) { #factor.level in factor.list) { factor.level <- factor.list[i] lines(pred.frame[pred.frame[,1]==factor.level,2], prediction[pred.frame[,1]==factor.level], lty = i) } # 现在绘制图例 if(is.null(leg.coords)){ x.max <- max(pred.frame[,2]) x.min <- min(pred.frame[,2]) x.range <- x.max - x.min x.pos <- c(x.min + (0.02 * x.range),x.min + (0.3 * x.range)) y.max <- max(prediction) y.min <- min(prediction) y.range <- y.max - y.min y.pos <- c(y.min + (0.8 * y.range),y.min + (0.95 * y.range)) legend(x = x.pos, y = y.pos, factor.list, lty = c(1:length(factor.list)), bty = "n") } else { legend(x = leg.coords[1], y = leg.coords[2], factor.list, lty = c(1:length(factor.list)), bty = "n", ncol = 2) } }}