完整数据集包含~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]'会导致代码运行。
我已经浏览了以下主题:
Possible directed cycle error in JAGS
它们中的问题似乎是由等式两边的“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的比例)。
否则你将不得不以某种方式重新制定你的模型以打破定向循环。这将涉及考虑系统底层的生物过程以及如何在不创建定向循环的情况下表示您想要的关系,因此我们无法提供帮助。
希望有所帮助,
马特