JAGS错误-可能涉及以下部分或所有节点的定向循环

Question

整个数据集包含大约11,000行。我一直在用K=400运行代码，同时检查代码是否运行。

所有行都与地图上的特定单元相关，并包含从Sentinel-2图像和数字高程地图中提取的信息。

117个细胞的子集还包含野外旅行中记录的栖息地协变量。因此，包括响应变量(S1和S2)和tussac的一些列的特征在于高比例的NAs。

代码：

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*S1i‘时，我得到以下错误：

Possible directed cycle involving some or all of the following nodes

然后列出S1、S2、lambda1和lambda2的所有值。

删除'c4*S1i‘会导致代码运行。

我已经看过以下的帖子了：

Possible directed cycle error in JAGS

https://stats.stackexchange.com/questions/220312/coding-a-jags-error-model-for-a-dependent-variable-that-has-increasing-variance

它们中的问题似乎是由于海报在等式的两边都使用了'y‘造成的。我认为我的问题是因为a4将代码发送到S2部分，然后c4将其发送回S1部分，这有点像定向循环。你知道怎么解决这个问题吗？

我已经包含了dataset的顶行，以防它有任何用处：

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

Answer 1

正如您已经正确识别的，您的问题是模型图中的一个有向循环。这是一个问题的原因是，事实证明，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)，否则此技巧将不起作用。

否则，您将不得不以某种方式重新制定您的模型，以打破定向循环。这将涉及到思考你的系统底层的生物过程，以及如何在不创建定向循环的情况下表示你想要的关系，所以这并不是我们能真正帮助的东西。

希望这能有所帮助，

哑光