cycle - JAGS 错误 - 可能涉及以下部分或全部节点的定向循环
问题描述
完整的数据集包含约 11,000 行。我一直在使用 K=400 运行代码,同时检查代码是否运行。
所有行都与地图上的特定单元格相关,并包含从 Sentinel-2 图像和数字高程图中提取的信息。
117 个单元的子集还包含在实地考察中记录的栖息地协变量。因此,包括响应变量(S1 和 S2)和 tussac 在内的一些列的特点是 NA 比例很高。
编码:
add_c4 <- "model{
for(i in 1:K) {
S1[i]~dpois(lambda1[i])
lambda1[i]<-exp(a0+a1*DEM_slope[i]+a2*DEM_elevation[i]+a3*tussac[i]+a4*S2[i])
S2[i]~dpois(lambda2[i])
lambda2[i]<-exp(c0+c1*DEM_slope[i]+c2*DEM_elevation[i]+c3*tussac[i]+c4*S1[i])
muLogit_tussac[i]<-b0 + sentinel1[i] + sentinel3[i] + sentinel7[i] + sentinel8[i] + sentinel9[i] + DEM_slope[i]
Logit_tussac[i]~dnorm(muLogit_tussac[i], tau)
logit(tussac[i])<-Logit_tussac[i]
}
# Priors
a0~dnorm(0, 10)
a1~dnorm(0, 10)
a2~dnorm(0, 10)
a3~dnorm(0, 10)
a4~dnorm(0, 10)
b0~dnorm(0, 10)
b1~dnorm(0, 10)
b2~dnorm(0, 10)
b3~dnorm(0, 10)
c0~dnorm(0, 10)
c1~dnorm(0, 10)
c2~dnorm(0, 10)
c3~dnorm(0, 10)
c4~dnorm(0, 10)
tau~dgamma(0.001, 0.001)
#data# S1, S2, K, sentinel1, sentinel3, sentinel7, sentinel8, sentinel9, DEM_slope, DEM_elevation
#inits# a0, a2, a3, a4, b0, b1, b2, b3, c0, c2, c3, c4
#monitor# a0, a1, a2, a3, a4, b0, b1, b2, b3, tau, ped, dic, c0, c1, c2, c3, c4
}"
当我包含 'c4*S1[i]' 时,我收到以下错误:
Possible directed cycle involving some or all of the following nodes
然后继续列出 S1、S2、lambda1 和 lambda2 的所有值。
删除 'c4*S1[i]' 会导致代码运行。
我浏览了以下线程:
其中的问题似乎是由于海报在等式两边都使用了“y”引起的。我认为我的问题是由于 a4 将代码发送到 S2 部分,而 c4 将其发送回 S1 部分,这有点像定向循环。知道如何解决这个问题吗?
我已经包含了数据集的顶部行,以防万一它有任何用途:
S1 S2 Logit_tussac moisture DEM_slope DEM_aspect DEM_elevation sentinel1 sentinel2 sentinel3 sentinel4 sentinel5 sentinel6 sentinel7 sentinel8 sentinel9 sentinel10
NA NA NA NA 2.434239 168.5011 0.588606366 0.0413 0.0499 0.0531 0.1035 0.1862 0.1968 0.1808 0.1318 0.0400 0.0199
NA NA NA NA 3.705001 178.1289 1.007037127 0.0966 0.1108 0.1212 0.0855 0.0917 0.1063 0.0937 0.1842 0.0341 0.0161
NA NA NA NA 5.006181 180.0000 1.883010797 0.1309 0.1472 0.1361 0.0855 0.0917 0.1063 0.0937 0.1572 0.0341 0.0161
NA NA NA NA 5.006181 180.0000 2.758984468 0.0542 0.0512 0.0472 0.0145 0.0127 0.0092 0.0166 0.0510 0.0148 0.0080
数据集子集,以便仅包含远程和本地感测数据的 117 行:
S1 S2 Logit_tussac moisture DEM_slope DEM_aspect DEM_elevation sentinel1 sentinel2 sentinel3 sentinel4 sentinel5 sentinel6 sentinel7 sentinel8 sentinel9 sentinel10
NA NA NA NA 14.917334 256.1612 12.24432 0.0513 0.0588 0.0541 0.1145 0.1676 0.1988 0.1977 0.1658 0.1566 0.0770
0 0 -9.210240 1 23.803741 225.1231 16.88028 0.1058 0.1370 0.2139 0.2387 0.2654 0.2933 0.3235 0.2928 0.3093 0.1601
NA NA NA NA 20.789165 306.0945 18.52480 0.0287 0.0279 0.0271 0.0276 0.0290 0.0321 0.0346 0.0452 0.0475 0.0219
NA NA -9.210240 1 6.689442 287.9641 36.08975 0.0462 0.0679 0.1274 0.1535 0.1797 0.2201 0.2982 0.2545 0.4170 0.2252
0 0 -9.210240 1 25.476444 203.0659 23.59964 0.0758 0.1041 0.1326 0.1571 0.2143 0.2486 0.2939 0.2536 0.3336 0.1937
1 0 -1.385919 3 1.672511 270.0000 39.55215 0.0466 0.0716 0.1227 0.1482 0.2215 0.2715 0.3334 0.2903 0.3577 0.1957
解决方案
正如您已经正确识别的那样,您的问题是模型图中的有向循环。这是一个问题的原因是,DAG(有向无环图)不包含任何有向环非常重要,否则无法保证我们可以定义稳定的后验来进行采样。
例如,采用以下包含有向循环的模型:
model <- 'model{
for(i in 1:N){
a[i] ~ dnorm(b[i], tau)
b[i] ~ dnorm(a[i], tau)
}
tau ~ dgamma(0.01,0.01)
#monitor# tau
#data# N
}'
N <- 10
runjags::run.jags(model)
没有明智的方法来估计这个模型,JAGS 会告诉你。但理论上可以估计这个模型:
model <- 'model{
for(i in 1:N){
a[i] ~ dnorm(b[i], tau)
b[i] ~ dnorm(a[i], tau)
}
tau ~ dgamma(0.01,0.01)
#monitor# tau
#data# N, a
}'
N <- 10
a <- rnorm(N)
runjags::run.jags(model)
改变的是所有 a[] 现在都是固定的(观察到的),所以我们实际上可以估计这个模型。但是 JAGS 仍然会检测到有向循环,因此需要一种解决方法:
model <- 'model{
for(i in 1:N){
a[i] ~ dnorm(b[i], tau)
b[i] ~ dnorm(aa[i], tau)
}
tau ~ dgamma(0.01,0.01)
#monitor# tau
#data# N, a, aa
}'
N <- 10
a <- rnorm(N)
aa <- a
runjags::run.jags(model)
这通过欺骗 JAGS 认为 a[] 和 aa[] 不相关来隐藏定向循环。但这仅在观察/修复所有 a[] 时才有效,否则模型中不会估计或定义丢失的 aa[]。在您的情况下,似乎部分观察到 S1[] 和 S2[],因此除非您可以简单地省略缺少 S1 或 S2 的行/观察值(这可能不可行,因为您说它们具有高NA的比例)。
否则,您将不得不以某种方式重新制定模型以打破定向循环。这将涉及考虑您系统背后的生物学过程以及如何在不创建定向循环的情况下表示您想要的关系,因此我们无法真正帮助解决。
希望有帮助,
马特
推荐阅读
- jenkins - 如果自上次提交以来某个目录中的任何文件已被修改(不是 FS 触发器),让 Jenkins 执行构建步骤
- swiftui - 如何为 SwiftUI 制作滑动手势?
- c# - 数据表计算时,对象不能从 DBNull 转换为其他类型
- python - 使用另一个列表将列表解压缩到单独的列表中
- gekko - Gekko的工业应用
- bash - 在 MacOS 上的 bash 脚本中使用关联数组
- r - 如何对包含数字的字符向量进行排序
- docker - Hazelcast 分布式缓存是否适用于 Docker Swarm
- mongoose - 猫鼬的连接是特定于版本的吗?
- java - 管理 ImageIcons 的 Java 最佳实践