从多个连续变量制作回归估计的曲面图

Question

我有一个包含分类和连续变量和样条的多层次模型。美好而复杂。无论如何，我正在尝试可视化模型拟合。

例如，这里有一些玩具数据：

library(lme4)
library(rms)
library(gridExtra)

## Make model using sleepstudy data
head(sleepstudy)
# Add some extra vars
sleepstudy$group <- factor( sample(c(1,2), nrow(sleepstudy), replace=TRUE) )
sleepstudy$x1 <- jitter(sleepstudy$Days, factor=5)^2 * jitter(sleepstudy$Reaction)

# Set up a mixed model with spline
fm1 <- lmer(Reaction ~ rcs(Days, 4) * group + (rcs(Days, 4) | Subject), sleepstudy)
# Now add continuous covar
fm2 <- lmer(Reaction ~ rcs(Days, 4) * group + x1 + (rcs(Days, 4) | Subject), sleepstudy)

# Plot fit
new.df <- sleepstudy
new.df$pred1 <- predict(fm1, new.df, allow.new.levels=TRUE, re.form=NA)
new.df$pred2 <- predict(fm2, new.df, allow.new.levels=TRUE, re.form=NA)

g1 <- ggplot(data=new.df, aes(x=Days)) +
    geom_line(aes(y=pred1, col=group), size=2) +
    ggtitle("Model 1")
g2 <- ggplot(data=new.df, aes(x=Days)) +
    geom_line(aes(y=pred2, col=group), size=2) +
    ggtitle("Model 2")

 grid.arrange(g1, g2, nrow=1)

绘图 1 是平滑的，但由于 x1 的影响，绘图 2 是锯齿状的。所以我想用 x = Days、y = x1 和 z = pred2 制作一个表面图，并按 group 分层。没有曲面图的经验，我开始使用线框命令：

wireframe(pred2 ~ Days * x1, data = new.df[new.df$group==1,],
      xlab = "Days", ylab = "x1", zlab="Predicted fit"
)

然而，虽然这个命令没有报错，但我的绘图是空白的：

问题：

我的线框哪里出错了？
有没有更好的方法来可视化我的模型拟合？

Answer 1

我发现 wireframe' orplot_ly' 表面所需的数据格式是 x 行乘以对应 z 值的 y 列的二维矩阵（我得到了对此的提示来自这个问题）。我还意识到我可以使用“expand.grid”来制作一个覆盖可能的 x 和 y 值范围的矩阵，并使用它们来预测 z，如下所示：

days <- 0:9
x1_range <- range(sleepstudy$x1)[2] * c(0.05, 0.1, 0.15, 0.2, 0.25, 0.3)

new.data2 <- expand.grid(Days = days, x1 = x1_range, group = unique(sleepstudy$group) )
new.data2$pred <- predict(fm2, new.data2, allow.new.levels=TRUE, re.form=NA)

然后我可以将它们填充到两个不同的矩阵中以表示模型中每个组的 z 面：

surf1 <- ( matrix(new.data2[new.data2$group == 1, ]$pred, nrow = length(days), ncol = length(x1_range)) )
surf2 <- ( matrix(new.data2[new.data2$group == 2, ]$pred, nrow = length(days), ncol = length(x1_range)) )
group <- c(rep(1, nrow(surf1)), rep(2, nrow(surf2) ))

最后我可以使用 plot_ly 绘制每个表面：

plot_ly (z=surf1, x = mets_range, y = ages, type="surface") %>%
    add_surface (z = surf2, surfacecolor=surf2,
             color=c('red','yellow'))

结果图：

所以结果图就是我想要的（虽然在这个虚构的例子中不是很有用，但在真实数据中很有用）。我唯一想不通的是如何显示两种不同的色标。我可以完全抑制比例，但如果有人知道如何为不同的表面显示 2 个比例，请告诉我，我将编辑答案。

从多个连续变量制作回归估计的曲面图

Making surface plot of regression estimates from multiple continuous variables

plot

model

r

surface

问题：