OpenBUGS 错误未定义变量

Question

我正在使用 OpenBUGS 和 R 包研究二项式混合模型R2OpenBUGS。我已经成功构建了更简单的模型，但是一旦我添加另一个级别来进行不完善的检测，我就会不断收到错误variable X is not defined in model or in data set。我尝试了许多不同的方法，包括更改数据结构以及将数据直接输入 OpenBUGS。我发布此内容是希望其他人有此错误的经验，并且也许知道为什么 OpenBUGS 无法识别变量 X，尽管据我所知它已明确定义。

我也收到了错误expected the collection operator c error pos 8- 这不是我以前收到的错误，但我也同样感到困惑。

模型和数据模拟函数均来自 Kery 的《生态学家 WinBUGS 简介》(2010)。我要注意的是，这里的数据集是代替我自己的数据的，这是相似的。

我包括构建数据集和模型的函数。抱歉长度。

# Simulate data: 200 sites, 3 sampling rounds, 3 factors of the level 'trt', 
# and continuous covariate 'X'

data.fn <- function(nsite = 180, nrep = 3, xmin = -1, xmax = 1, alpha.vec = c(0.01,0.2,0.4,1.1,0.01,0.2), beta0 = 1, beta1 = -1, ntrt = 3){
  y <- array(dim = c(nsite, nrep))  # Array for counts
  X <- sort(runif(n = nsite, min = xmin, max = xmax))   # covariate values, sorted
  # Relationship expected abundance - covariate
  x2 <- rep(1:ntrt, rep(60, ntrt)) # Indicator for population
  trt <- factor(x2, labels = c("CT", "CM", "CC"))
  Xmat <- model.matrix(~ trt*X)
  lin.pred <- Xmat[,] %*% alpha.vec # Value of lin.predictor
  lam <- exp(lin.pred)
  # Add Poisson noise: draw N from Poisson(lambda)
  N <- rpois(n = nsite, lambda = lam)
  table(N)                # Distribution of abundances across sites
  sum(N > 0) / nsite          # Empirical occupancy
  totalN <- sum(N)  ;  totalN
  # Observation process
  # Relationship detection prob - covariate
  p <- plogis(beta0 + beta1 * X)
  # Make a 'census' (i.e., go out and count things)
  for (i in 1:nrep){
    y[,i] <- rbinom(n = nsite, size = N, prob = p)
  }
  # Return stuff
  return(list(nsite = nsite, nrep = nrep, ntrt = ntrt, X = X, alpha.vec = alpha.vec, beta0 = beta0, beta1 = beta1, lam = lam, N = N, totalN = totalN, p = p, y = y, trt = trt))
}

data <- data.fn()

这是模型：

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
    y[i,j] ~ dbin(p[i,j], C[i])
    lp[i,j] <- beta0 + beta1*X[i]
    p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(C = y, trt = as.numeric(trt), X = s.X)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

Answer 1

注意：在我注意到代码的另一个问题之后，这个答案已经进行了重大修改。

如果我正确理解你的模型，那么你会混淆模拟数据中的y和N ，以及作为C传递给 Bugs 的内容。您将y变量（矩阵）传递给 Bugs 模型中的 C 变量，但这是作为向量访问的。据我所知，C代表二项式绘制（实际丰度）中的“试验”次数，即数据集中的N 。变量y（矩阵）在模拟数据和 Bugs 模型中被称为相同的东西。

据我了解，这是对您的模型的重新表述，并且运行正常：

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
        y[i,j] ~ dbin(p[i,j], C[i])
        lp[i,j] <- beta0 + beta1*X[i]
        p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
N<- data$N
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(y = y, trt = as.numeric(trt), X = s.X, C= N)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

总体而言，该模型的结果看起来不错，但 beta0 和 beta1 存在较长的自相关滞后。beta1 的估计似乎也有点偏差（~= -0.4），因此您可能需要重新检查 Bugs 模型规范，以便它与模拟模型匹配（即您正在拟合正确的统计模型）。目前，我不确定是否如此，但我现在没有时间进一步检查。