我正在改编一些开源代码来研究环境和有机变量如何影响感染概率。我收到一个错误，指出第 10 行可能存在涉及 resid_pred 的定向循环。很抱歉代码很长，但我想包含所有内容，因为错误可能发生在模型中的其他地方？非常感谢任何帮助。
cat("
model {
#Sensitivity and specificity
spec <- sum(ifelse(equals(pred, infected) && infected == 0, 1, 0)) / sum(equals(infected, 0))
sens <- sum(ifelse(equals(pred, infected) && infected == 1, 1, 0)) / sum(equals(infected, 1))
#Posterior predictive check
ppc <- step(fit_pred - fit_data)
fit_data <- sum(resid_data)
fit_pred <- sum(resid_pred)
resid_data <- abs(infected) - prob
resid_pred <- abs(pred) - prob
#Phylogenetic logistic regression model
#Logistic regression w/ varying intercepts for each host species
for(i in 1:n_obs){
#likelihood
infected[i] ~ dbern(prob[i])
#simulated prediction
pred[i] ~ dbern(prob[i])
#linear predictor
logit(prob[i]) <- reg[ecoregion[i]] + alpha[hostsp[i]] +
phy[hostsp[i]] + beta_male * male[i] +
beta_bodymass * bodymass[i] + beta_migratory * migratory[i] + beta_adult * adult[i] +
beta_high_canopy * high_canopy[i] + beta_min_rain * min_rain[i] +
beta_max_rain * max_rain[i] + beta_min_temp * min_temp[i] +
beta_max_temp * max_temp[i] + beta_ndvi * ndvi[i] +
beta_ndwi * ndwi[i] +
beta_elevation * elevation[i]
}
#Impute any missing covariate values
for(i in 1:n_obs){
male[i] ~ dbern(0.5);
adult[i] ~ dbern(0.5)
}
#Estimate varying intercepts among regions
for(r in 1:n_ecoregion){
reg[r] ~ dnorm(0, tau_reg);
}
tau_reg <- pow(sigma_reg,-2); sigma_reg ~ dexp(0.5)
#Estimate varying intercepts among host species from normal distributions
for (j in 1:n_species){
alpha[j] ~ dnorm(0, tau_alpha);
}
tau_alpha <- pow(sigma_alpha,-2); sigma_alpha ~ dexp(0.5)
#Multivariate normal phylogenetic component
#invA_array = array of inverse phylo vcv matrices
#K indicates which vcv matrix to sample (see below)
phy[1:n_species] ~ dmnorm(zeros[], tau_phy*invA_array[,,K])
#Precision for phylogenetic cov matrix (inverted bc phylo vcv matrices are inverted)
tau_phy <- 1 / pow(sigma_phy, 2); sigma_phy ~ dexp(0.5)
#Prior for sampling phylogenetic vcv matrices with equal probabilities
K ~ dcat(p[])
for (k in 1:n_tree) {
p[k] <- 1 / n_tree
}
#Gibbs variable selection for beta coefficients
beta_high_canopy <- beta_high_canopyT * Ind_high_canopy
beta_high_canopyT ~ dnorm(0, tau_sel[1, Ind_high_canopy + 1])
Ind_high_canopy ~ dbern(pind) #indicators determine whether to sample from spike or slab
beta_min_rain <- beta_min_rainT * Ind_min_rain
beta_min_rainT ~ dnorm(0, tau_sel[2, Ind_min_rain + 1])
Ind_min_rain ~ dbern(pind)
beta_max_rain <- beta_max_rainT * Ind_max_rain
beta_max_rainT ~ dnorm(0, tau_sel[3, Ind_max_rain + 1])
Ind_max_rain ~ dbern(pind)
beta_min_temp <- beta_min_tempT * Ind_min_temp
beta_min_tempT ~ dnorm(0, tau_sel[4, Ind_min_temp + 1])
Ind_min_temp ~ dbern(pind)
beta_max_temp <- beta_max_tempT * Ind_max_temp
beta_max_tempT ~ dnorm(0, tau_sel[5, Ind_max_temp + 1])
Ind_max_temp ~ dbern(pind)
beta_male <- beta_maleT * Ind_male
beta_maleT ~ dnorm(0, tau_sel[6, Ind_male + 1])
Ind_male ~ dbern(pind)
beta_bodymass <- beta_bodymassT * Ind_bodymass
beta_bodymassT ~ dnorm(0, tau_sel[7, Ind_bodymass + 1])
Ind_bodymass ~ dbern(pind)
beta_ndvi <- beta_ndviT * Ind_ndvi
beta_ndviT ~ dnorm(0, tau_sel[8, Ind_ndvi + 1])
Ind_ndvi ~ dbern(pind)
beta_ndwi<-beta_ndwiT*Ind_ndwi
beta_ndwiT~dnorm(0,tau_sel[9,Ind_ndwi+1])
Ind_ndwi~dbern(pind)
beta_migratory<-beta_migratoryT*Ind_migratory
beta_migratoryT~dnorm(0,tau_sel[10,Ind_migratory+1])
Ind_migratory~dbern(pind)
beta_adult<-beta_adultT*Ind_adult
beta_adultT~dnorm(0,tau_sel[11,Ind_adult +1])
Ind_adult ~dbern(pind)
beta_elevation<-beta_elevationT*Ind_elevation
beta_elevationT~dnorm(0,tau_sel[12,Ind_elevation+1])
Ind_elevation~dbern(pind)
#Sample the indicators and formulate the spike and slab distributions
for(s in 1:12){
tau_sel[s,2]<-pow(sd_sel[s],-2);
#slab is a vague normal prior
sd_sel[s]~dexp(0.5);
tau_sel[s,1]~dunif(85,100);
#spike concentrates at zero with low variance
}
pind <-0.2
#coefficients have ~ 20% inclusion probability
}
",
file=(phylo.jags.mod<-tempfile()))
#Model Run Parameters
na = 25000 # adaptation
nb = 3e+05 # burnin
ni = 25000 # samples
nc = 2 # chains
library(MCMCpack)
library(runjags)
Sys.setenv(JAGS_HOME = "C:/Program Files/JAGS/JAGS-4.3.1")
library(rjags)
load.module("glm")
param<-c("reg","alpha",
"beta_elevation","beta_ndwi","beta_ndvi",
"beta_high_canopy","beta_male",
"beta_migratory","beta_min_rain",
"beta_max_rain","beta_min_temp",
"beta_adult","beta_max_temp",
"beta_bodymass","Ind_elevation",
"Ind_adult","Ind_high_canopy",
"Ind_male","Ind_migratory","Ind_min_rain",
"Ind_bodymass","Ind_max_rain",
"Ind_min_temp","Ind_max_temp",
"Ind_ndwi","Ind_ndvi","phy",
"sigma_reg","sigma_alpha","sigma_phy",
"ppc","spec","sens")
jags_inits_haem <- function() {
list(reg = rnorm(data_jags_haem$n_ecoregion),
alpha = rnorm(data_jags_haem$n_species),
phy = rnorm(data_jags_haem$n_species),
beta_elevationT = rnorm(1),
beta_ndwiT = rnorm(1),
beta_ndviT = rnorm(1),
beta_high_canopyT = rnorm(1),
beta_min_rainT = rnorm(1),
beta_max_rainT = rnorm(1),
beta_maleT = rnorm(1),
beta_migratoryT = rnorm(1),
beta_adultT = rnorm(1),
beta_min_tempT = rnorm(1),
beta_max_tempT = rnorm(1),
beta_bodymassT = rnorm(1),
sigma_reg = runif(1, 0.01,10),
sigma_alpha = runif(1,0.01, 10),
sigma_phy = runif(1,0.01, 10),
sd_sel = runif(12,0.01, 10))
}
load("FOLDER/phyloglmm_data_jags_haem.rda")
rjags_phylo_haem <- jags.model(phylo.jags.mod,
data = data_jags_haem, inits = jags_inits_haem,
n.chains = nc, n.adapt = na)
update(rjags_phylo_haem, nb)
out_haem <- jags.samples(rjags_phylo_haem,
param, ni, thin = 25)
save(out_haem, file = "Processeddata/out_haem.rda")
我读过其他帖子，但我无法确定定向循环来自哪里。我也在 SourceForge 中发布了此内容，其中包含“phyloglmm_data_jags_haem.rda”文件的附件。我还没弄清楚如何将其附加到这里 https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/9c6dc6b0ef/
我可以运行你的代码。 我对结果的准确性不做任何承诺。
问题似乎是您在输入数据的第 326 次观察中缺少预测变量
ndvifor(i in 1:n_obs){
male[i] ~ dbern(0.5);
adult[i] ~ dbern(0.5)
ndvi[i] ~ dunif(1100, 8000)
}
允许模型编译。 [1100 和 8000 分别略低于/高于
ndvi顺便说一句，为了打发我的无聊，我将你的 MCMC 参数调整为
na = 250 # adaptation
nb = 3e+02 # burnin
ni = 250 # samples
nc = 1 # chains
这允许进行健全性检查，但对于实际生产工作来说显然太小了。
我最好的猜测是，缺失的
ndvindviprob顺便说一句，您不需要
library(MCMCpack)我是如何解决这个问题的？
我首先注释掉似然循环内的所有内容。 这表明先验是可估计的。 然后，我又取消了预测模型中每个术语的注释。 每个模型都已编译且可估算，直到我将其包含在内
beta_ndvi * ndvi[i]在做其他事情之前检查所有预测变量中是否有缺失值会更快！