Ошибка 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]' приводит к выполнению кода.
Я просмотрел следующие темы:
Возможная направленная ошибка цикла в 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
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).
В противном случае вам придется как-то переформулировать свою модель, чтобы разорвать направленный цикл. Это будет включать в себя размышления о биологическом процессе, лежащем в основе вашей системы, и о том, как представлять отношения, которые вы хотите, без создания направленного цикла, так что на самом деле мы не можем помочь с этим.
Надеюсь, это поможет,
Matt