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

Question

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

所有行都与 map 上的特定单元格相关，并包含从 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

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

其中的问题似乎是由张贴者在等式两边都使用“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

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 的比例)。

否则您将不得不以某种方式重新制定您的模型以打破定向循环。这将涉及考虑系统背后的生物过程，以及如何在不创建定向循环的情况下表示您想要的关系，因此这不是我们真正可以提供帮助的事情。

希望对你有帮助，

马特