r - OpenBUGS 错误 undefined variable

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 的 Introduction to WinBUGS for Ecoologists (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 传递给 Bug 的内容混淆了。您将 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 模型规范，以便它与模拟模型匹配(即您正在拟合正确的统计模型)。目前，我不确定它是否确实如此，但我现在没有时间进一步检查。