4d 曲面图逻辑回归（x、y、z、颜色）

Question

我正在尝试绘制涉及 4 个定量变量（3 个预测变量和 1 个响应）的二项式关系的 4d 曲面图。 我有森林覆盖率 (cov)、景观异质性 (het) 和大量捕食者 (abund) 来解释捕食的概率 (pred)。 所以我想在图表上放 x = cov, y = het, z = abund 和一个颜色渐变来代表pred。我想要一个类似这样的图表：

我的模型是：

model <- glm(cbind(pred,npred)~ cov + het + abund + abund:cov + cov:het, data=results, family=binomial)

我绘制图表的脚本：

het <- seq(0, 2, 0.05288)
cov <- seq(0, 0.8, 0.021600)
abun <- seq(0, 0.9, 0.024)

gg <- expand.grid(het=het, cov=cov, abund=abund)
gg$pred < -predict(model,newdata=gg, type = "response")
head(gg)

with(gg, plot3d(het, cov, abund, type = "n"))
cols <- heat.colors(20)
cuts <- with(gg, cut(pred, breaks = 20))

with(gg, plot3d(het, cov, abund, type = "n"))
with(gg, surface3d(1:20,1:20, matrix(pred, ncol = 20),
                  color = cols[cuts], back = "fill"))

不幸的是，我认为有些地方出了问题，并没有像我想象的那样情节。 有人可以帮我编写脚本或提出更好的建议吗？

我的数据集：

het cov abun    pred    npred
1.90164 0.636460117 0.730424235 24  14
2.13681 0.409547496 0.735686442 8   42
2.18218 0.631526938 0.744585387 2   48
2.10223 0.568536934 0.746072471 8   36
2.11588 0.390850929 0.747720518 2   50
2.10288 0.465843118 0.754245798 6   42
2.28891 0.272948757 0.757754057 4   50
2.31959 0.281255571 0.776337246 2   32
2.10626 0.426842961 0.778512783 4   46
2.0163  0.546051041 0.780669728 6   50
1.91282 0.433635614 0.809807372 2   48
1.67011 0.594238535 0.827360862 18  26
1.80076 0.719426666 0.830916931 4   26
2.29307 0.471266098 0.839754622 2   44
2.00811 0.574829217 0.839935299 2   40
2.25667 0.416338714 0.849465507 0   54
2.02098 0.306265093 0.860655587 4   42
2.10236 0.263078082 0.860656944 6   36
1.23721 0.236391255 0.879960428 6   44

Answer 1

您的预测是三个预测变量的函数，因此您不能将它们表示为单个表面。有两种常见的方法可以做到这一点：

1. 使用 misc3d::contour3d 函数绘制 3d 轮廓。

2. 为变量之一的不同值绘制多个曲面。

在执行第一个操作之前，需要修复代码中的一些拼写错误：有时您使用abun，有时使用abund：我将始终使用abun。你也有

gg$pred < -predict(model,newdata=gg, type = "response")

<- 中间的空格完全改变了含义。最后，使用dput(results) 可以让其他人更容易使用数据。一切都清理干净了：

results <- structure(list(het = c(1.90164, 2.13681, 2.18218, 2.10223, 2.11588, 
2.10288, 2.28891, 2.31959, 2.10626, 2.0163, 1.91282, 1.67011, 
1.80076, 2.29307, 2.00811, 2.25667, 2.02098, 2.10236, 1.23721
), cov = c(0.636460117, 0.409547496, 0.631526938, 0.568536934, 
0.390850929, 0.465843118, 0.272948757, 0.281255571, 0.426842961, 
0.546051041, 0.433635614, 0.594238535, 0.719426666, 0.471266098, 
0.574829217, 0.416338714, 0.306265093, 0.263078082, 0.236391255
), abun = c(0.730424235, 0.735686442, 0.744585387, 0.746072471, 
0.747720518, 0.754245798, 0.757754057, 0.776337246, 0.778512783, 
0.780669728, 0.809807372, 0.827360862, 0.830916931, 0.839754622, 
0.839935299, 0.849465507, 0.860655587, 0.860656944, 0.879960428
), pred = c(24L, 8L, 2L, 8L, 2L, 6L, 4L, 2L, 4L, 6L, 2L, 18L, 
4L, 2L, 2L, 0L, 4L, 6L, 6L), npred = c(14L, 42L, 48L, 36L, 50L, 
42L, 50L, 32L, 46L, 50L, 48L, 26L, 26L, 44L, 40L, 54L, 42L, 36L, 
44L)), class = "data.frame", row.names = c(NA, -19L))

model <- glm(cbind(pred,npred)~ cov + het + abun + abun:cov + cov:het,
             data=results, family=binomial)

het <- seq(0, 2, 0.05288)
cov <- seq(0, 0.8, 0.021600)
abun <- seq(0, 0.9, 0.024)
gg <- expand.grid(het=het, cov=cov, abun=abun)
gg$pred <- predict(model,newdata=gg, type = "response")
head(gg)
#       het cov abun         pred
# 1 0.00000   0    0 2.220446e-16
# 2 0.05288   0    0 2.220446e-16
# 3 0.10576   0    0 2.220446e-16
# 4 0.15864   0    0 2.220446e-16
# 5 0.21152   0    0 2.220446e-16
# 6 0.26440   0    0 2.220446e-16

要绘制等高线图，您需要选择等高线水平。您可以使用

查看预测范围

range(gg$pred)
# [1] 2.220446e-16 1.000000e+00

所以你的模型预测的概率基本上是 0 到 1。我想要 0 到 1 之间的 20 个等高线，但 contour3d 函数不喜欢超出范围 等高线，所以我将使用 1/20，...，19/20：

levels <- (1:19)/20
colors <- heat.colors(19)

contour3d 函数需要数组中的值，而不是向量：

pred <- array(gg$pred, c(length(het), length(cov), length(abun)))

现在画一些轴，然后是情节：

plot3d(gg, type = "n")
contour3d(pred, levels, het, cov, abun, color = colors, alpha = 0.5, add = TRUE)

alpha = 0.2 值使表面部分透明；您可能想尝试 0 到 1 之间的其他值，看看您喜欢什么。这就是我得到的。在 R 中，您可以旋转绘图以从不同角度查看它，或调用 rglwidget() 以将其保存为可在网络浏览器上查看的可旋转形式。

我会留给你解决如何做第二种情节。