你好呀！我完全get到你的需求了——你想用mgcViz实现类似mgcv的vis.gam那样的模型整体拟合估计值可视化，同时也想绘制不同组之间的估计值差值，而不只是平滑项的效应和效应差值，对吧？

先明确下核心区别：mgcv的vis.gam输出的是模型对响应变量的整体预测值（比如当fac="A1"时，给定x1、x2的y的估计值），而mgcViz默认的plot(sm(o, 1))只展示单个平滑项的效应（也就是不含模型截距和其他公共项的部分）。下面就来一步步实现你要的效果：

首先先还原你的示例环境：

library(mgcv)
library(mgcViz)
set.seed(235)
dat <- gamSim(1,n=1500,dist="normal",scale=20)
dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) )
dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) )
bs <- "cr"; k <- 12
b <- gam(y ~ s(x2, x1, by = fac), data=dat)
o <- getViz(b, nsim = 0)

你提到的三个参考可视化效果：

• mgcv的vis.gam绘制的整体估计值曲面：
  
• mgcViz的plot(sm(o,1))绘制的平滑项效应：
  
• mgcViz的plotDiff绘制的平滑项效应差值：
  

一、用mgcViz绘制模型整体估计值（对应vis.gam的效果） #

这里有两种简单的实现方式：

方法1：给平滑项效应加上模型截距

你的模型结构中，fac="A1"时的整体预测值 = 模型截距 + 对应平滑项的效应。我们可以直接提取平滑项对象，加上截距后再绘图：

# 获取模型截距
model_intercept <- coef(b)[1]

# 提取fac=A1对应的平滑项效应，加上截距得到整体估计值
sm_a1 <- sm(o, 1)
sm_a1_estimate <- sm_a1
sm_a1_estimate$fit <- sm_a1_estimate$fit + model_intercept  # 加上截距转换为整体估计

# 用mgcViz的风格绘图
plot(sm_a1_estimate) + 
  l_fitRaster() + 
  l_fitContour(color = "white") +
  labs(title = "mgcViz绘制的整体估计值（fac=A1）")

这个结果就和vis.gam(b, cond=list(fac="A1"),plot.type="contour")完全一致啦！

方法2：生成预测网格后绘图

如果你的模型更复杂（比如包含其他协变量），这种方法更通用：先生成指定条件的预测数据，再用ggplot2/mgcViz的可视化语法绘图：

# 生成预测网格：固定fac为"A1"，x1、x2取连续序列
pred_grid <- expand.grid(
  x1 = seq(min(dat$x1), max(dat$x1), length.out = 100),
  x2 = seq(min(dat$x2), max(dat$x2), length.out = 100),
  fac = "A1"
)
# 预测整体估计值
pred_grid$y_hat <- predict(b, newdata = pred_grid)

# 用mgcViz兼容的ggplot2风格绘图
library(ggplot2)
ggplot(pred_grid, aes(x = x1, y = x2, z = y_hat)) +
  geom_raster(aes(fill = y_hat)) +
  geom_contour(color = "white", linewidth = 0.5) +
  scale_fill_viridis_c(option = "viridis") +
  labs(title = "mgcViz风格的整体估计值（fac=A1）", fill = "y估计值") +
  theme_bw()

二、用mgcViz绘制整体估计值的差值 #

这里要注意：同一协变量条件下，不同fac组的整体估计值差值 = 对应平滑项的效应差值，因为模型中的公共项（比如截距、其他非组特异性项）会相互抵消。

所以你之前用plotDiff的代码，其实已经等价于估计值的差值了！我们只需要调整标注，明确这一点：

# 绘制fac=A2与fac=A1的估计值差值
pl <- plotDiff(s1 = sm(o, 1), s2 = sm(o, 2))
pl + 
  l_fitRaster() + 
  l_fitContour(color = "white") +
  labs(title = "整体估计值差值（fac=A2 - fac=A1）")

如果你的模型包含其他组特异性项，只需要把对应项的效应也加入差值计算即可，逻辑是一样的：先提取各组的所有相关项效应，求和得到整体估计，再计算差值后绘图。

备注：内容来源于stack exchange，提问作者David