r - 4d 曲面图逻辑回归(x、y、z、颜色)
问题描述
我正在尝试绘制涉及 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
解决方案
您的预测是三个预测变量的函数,因此您不能将它们表示为单个表面。有两种常见的方法可以做到这一点:
misc3d::contour3d
使用该函数 绘制 3d 轮廓。为变量之一的不同值绘制多个曲面。
在执行第一个操作之前,需要修复代码中的一些拼写错误:有时您使用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。我想要 20 个介于 0 和 1 之间的轮廓,但该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()
以将其保存为可在 Web 浏览器上查看的可旋转形式。
我会留给你来解决如何做第二种情节。
推荐阅读
- powershell - Switch 参数和 If 语句
- python-3.x - 使用 pandas python3 从特定行更新特定值
- html - 我不断收到“加载资源失败:net::ERR_FILE_NOT_FOUND”错误,不知道如何解决
- python - 字典列表:按内部字典键分组聚合值
- c++ - 阅读最后一行代码算法的问题 - C++
- javascript - Nuxt js 外部 URL 和自动关闭
- c++ - 有哪些有效的方法可以对一组 > 100 万个字符串进行重复数据删除?
- arrays - 尝试读取我的 firebase 数据库并创建包含其内容的对象数组时,我不断收到 SIGABRT 错误。我究竟做错了什么?
- php - 指令“track_errors”已弃用
- c++ - 以线程安全的方式调用 `QQuickPaintedItem::updateImage(const QImage& image)`(无 QThread)